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SUMMARY 


An improved computer program for studying the linear ocean wave refraction 
process has been developed. The program features random-access modular bathym- 
etry data storage to minimize computer resource requirements. The user can 
select from three bottom topography approximation techniques which provide 
varying degrees of bathymetry data smoothing. This allows the user to assess 
the sensitivity of computed results to the degree of smoothing chosen. Wave- 
ray patterns can be generated and plotted automatically with either a specified 
uniform initial (deepwater) ray spacing or a constrained final (nearshore) ray 
spacing. As an option, wave-crest patterns can be constructed and plotted by 
using a cubic polynomial to approximate individual crest segments between 
adjacent rays. 


INTRODUCTION 

The prime objective of the NASA Oceanology Program is development and 
demonstration of the feasibility of satellite systems such as SEASAT-A for the 
acquisition of global data on ocean conditions. The broad-scale data provided 
by such systems can be used to initialize analytical models which would then 
provide short-term, higher resolution forecasts of ocean conditions for ship 
routing and offshore activities. One such analytical tool is the wave refrac- 
tion model, which can predict changes in the speed and direction of waves as 
they travel from deep water across the continental shelf toward the coastline. 
Such changing patterns affect the distribution of wave energy (wave height) 
over the area through which the waves are traveling and along the shoreline. 

Early models, such as that developed in reference 1, used graphical tech- 
niques for constructing refraction diagrams. A wave refraction diagram is a 
map showing either a group of wave crests at a given time or the successive 
positions of a particular wave crest as it moves shoreward. The first models 
required the construction of each wave crest as it advanced toward shore, but 
later models required only the construction of wave rays which are locally per- 
pendicular to each crest. The wave-ray technique eliminated the time-consuming 
process of constructing each crest, at the expense of reducing the visual famil 
iarity of the end product. This graphical ray technique was later incorporated 
into a computer program by Griswold in 1963 (ref. 2), but the program's useful- 
ness was limited because a spatial grid of wave speeds was required in calcu- 
lating wave-ray paths. The refraction model developed by Wilson (ref. 3) used 
a grid of water depths which eliminated many of the problems associated with 
the wave-speed grid; however, the model could not provide wave height. The 
model developed by Dobson (ref. 4) succeeded in coupling wave-ray construction 
using a grid of water depths with a finite-difference solution to the wave 
intensity equation developed by Munk and Arthur (ref. 5) to obtain wave heights 
The Dobson model formed the basis for the model used at the Langley Research 
Center (LaRC) in a recent broad-scale analytical study of wave refraction in 
the Virginian Sea region of the mid-Atlantic continental shelf (ref. 6). 


Several modifications have been incorporated into the LaRC refraction 
model since the Virginian Sea study. Bathymetry data are stored and used in 
modular form by using random-access techniques to pass data modules as required 
from a peripheral mass-storage file to the central memory of the computer 
(ref. 7). This technique allows very large geographical regions to be studied 
with a minimum of computer resources but does not restrict the applicability 
of the model to smaller areas. Ocean bottom topography can be approximated by 
using either the quadratic least squares, the cubic least squares, or the con- 
strained bicubic interpolation approaches (ref. 8). With these three approxi- 
mation options, the user can select different degrees of input data smoothing 
and assess the sensitivity of computed results to the degree of smoothing 
selected. In the normal mode, the spatial position of the initial point on 
each wave ray in a particular set is computed automatically to achieve a uni- 
form deepwater ray density. In the previous version of the model, the spatial 
coordinates of each initial point were input on individual computer cards. By 
selecting the uniform deepwater ray density, the user can readily identify 
areas of relatively high or low energy on the resultant refraction diagram as 
areas of convergent or divergent ray groups. As an alternative, the user can 
select a control option which provides for automatic computation of deepwater 
ray spacing to satisfy specified constraints on nearshore ray density. Selec- 
tion of this option could provide increased directional information in the 
nearshore zone while sacrificing the easy identification of high or low energy 
regions. As a third option, wave crests can be constructed by using a cubic 
polynomial to approximate the crest segment between corresponding points on 
adjacent rays. Selection of this crest output option also results in selection 
of the aforementioned controlled nearshore ray density option. 

This paper presents a description of the modified LaRC wave refraction 
model as developed for use on the Control Data Cyber 173 or 175 computer under 
Network Operating System (NOS) 1.2. A central memory field length of 660008 
is required, along with the peripheral mass-storage file used for input of 
bathymetry data. A description, flow diagram, and source code listing are pro- 
vided for the main program and each subroutine. Program input and output 
parameters are listed and described. Sample graphic outputs are provided. A 
list of system subroutines used, a detailed description of the wave-crest 
approximation technique and the accompanying subroutine for controlling near- 
shore ray density, and a description of a utility program which can be used for 
modularizing a regional bathymetry array are included in the appendixes. 


SYMBOLS 

a-|,a 2 ,a 3 coefficients of cubic polynomial used to approximate wave-crest 
segment 

b ray separation distance, meters 

c wave speed, meters/second 

d 
g 
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water depth, meters 

acceleration due to gravity, meters/second^ 


H wave height, meters 

Kj, refraction coefficient 

Kg shoaling coefficient 

L wavelength, meters 

p coefficient in equation governing ray separation, defined by 

equation (10), seconds"^ 

q coefficient in equation governing ray separation, defined by 

equation (11), seconds"^ 

s spatial coordinate along wave ray 

t time, seconds 

X alongshore coordinate, nautical miles 

Xo,Xjji computational boundaries along x-axis 

y offshore coordinate, nautical miles 

yQ,ym computational boundaries along y-axis 

a ray angle with respect to x-axis, degrees 

3 ray separation factor 

ic ray curvature, radians/meter 

a radian wave frequency, radians/second 

0 angle between adjoining first-order wave-crest segments, degrees 
Subscripts: 

1 ray number index 

j index of computational point along ray 

o initial 

A prime denotes differentiation with respect to water depth d. 

PROBLEM DESCRIPTION 

The problem considered in the present paper is that of computing the 
refraction characteristics of ocean waves propagating from deep water across 
the continental shelf toward the coastline. The results of linear wave theory. 
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which assumes a progressive sinusoidal wave of constant frequency and small 
steepness, have been adopted; and steady-state conditions of constant tide 
height and no winds are assumed. The governing equations have been given in 
references 4 and 5 and are discussed briefly in this section. Also discussed 
is the method of solution, which is based primarily on the method presented in 
reference 4. Two additional features enhance the solution in terms of bathym- 
etry data storage requirements and alternate numerical techniques for approxi- 
mating ocean bottom topography. 

any wave ray (orthogonal to a wave crest) can be represented 


cos a ( 1 ) 


The path of 
by the equations 

dx 

— = c 
dt 


dy 

— = c sin a 
dt 


( 2 ) 


da 1 /3d 3d 

K = — = — c' — sin a - — 
ds c \3x 3y 


where d is the local water depth. Ray orientation relative to regional 
coordinates is depicted in the following sketch; 


cos a 


(3) 


y 



The local wave phase speed c can be expressed according to linear wave 
theory by 


ad 

c - Cq tanh — = 0 
c 


( 4 ) 
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where Cq is the initial (deepwater) phase speed and a is the fixed radian 
wave frequency. On the basis of properties of the hyperbolic tangent function, 
the quantity tanh ad/c in equation (4) can be replaced by 1 if the value 
of ad/c is sufficiently large (deepwater range), or by the quantity od/c 
if the value of ad/c is sufficiently small (shallow-water range). For the 
purposes of the present paper, deep water is assumed if the local water depth 
exceeds 0.25 times the initial wavelength (ref. 6). Shallow water is assumed 
if the local depth is less than 0.005 times the wavelength. All other depths 
lie within the intermediate range in which the full transcendental form of 
equation (4) must be used. By differentiating equation (4) with respect to d, 
it can be shown that 



(5) 


By assuming that energy is neither transmitted across wave rays nor dissipated 
by bottom friction or percolation, the wave height at any point along a wave 
ray can be computed by 


H = HoKgKp 


( 6 ) 


where Hq is the initial wave height. The local shoaling coefficient Kg can 
be expressed as (ref. 4) 


Kg = 


/ , _£a_d./c \ 

y sinh 2ad/c/ 


1/2 


(7) 


The local refraction coefficient Kp is given by 


Kr 



(B)-1/2 


( 8 ) 


The ray separation factor 3 is governed by the differential equation 


d23 d3 

— + p — + q3 = 0 
dt2 dt 


(9) 
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in which 


/9d 3d ' 

p = -2c ' — cos a + — sin a 
\3x 3y 


( 10 ) 


and it can be shown that 


q 


cc 



/ 32d 

- 2 sin a cos ot 

\3x 8y 


3d 9d\ 



( 11 ) 


where 



( 12 ) 


The variables of primary influence in the solution of the governing equa- 
tions are the local water depth and its spatial variation. For a large geo- 
graphical area such as the continental shelf, a practical method for the input 
of bathymetry data to the refraction model is through a spatial grid format in 
which known depth values are supplied at the nodes of the grid. Such a format 
is used in the present model, but instead of storing the bathymetry data array 
for the entire region of interest in central memory, as with the model described 
in reference 4, overlapping subregional modules of bathymetry data are read as 
required into memory from a peripheral mass-storage file by using nonsequential 
random-access techniques. This method of data storage is described in detail 
in reference 7 and was shown to decrease computer storage requirements by 
75 percent for a refraction model of the mid-Atlantic continental shelf. 

In conjunction with the grid- type format for bathymetry data, a technique 
is required for determining depth values at intermediate points between nodes 
and for approximating spatial derivatives of the depth as required in equa- 
tions (3), (10), and (11). The present model offers the selection of three 
such bottom topography approximation techniques, which are discussed in refer- 
ence 8. These include the quadratic least squares technique, which was used 
in the model described in reference 4, the cubic least squares technique, and 
the constrained bicubic interpolation technique. Both least squares techniques 
offer some smoothing of the input bathymetry data but do not guarantee conti- 
nuity in the computed depth or its derivatives. The constrained bicubic inter- 
polation technique guarantees continuity in the depth values and first-order 
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and mixed second-order partial derivatives; it has the possible disadvantage 
of using actual input data without smoothing. 

Apart from the techniques for bathymetry data storage and bottom topog- 
raphy approximation, the basic solution technique used in the present model is 
described in detail in reference The path of each wave ray is constructed 
in equal time steps by an iterative solution of equations (1) to (3). At each 
step the wave height is calculated by applying Fox's finite-difference solution 
(see ref. 4 for details) to equation (9) to determine the next ray separation 
factor and then using this result in equation (6). 


PROGRAM DESCRIPTION 

The LaRC wave refraction computer program is coded in FORTRAN Extended 
(FTN) Version 4 for use on the Control Data Cyber 173 or 175 computer. It con- 
sists of a main program WAVE and subroutines RAYINIT, RAYCON, REFRAC, CURVE, 
DEPTH, HEIGHT, WRITER, PLOTR, RAYOPT, and MATRIX, whose names are descriptive 
of the functions they perform. The functions performed by routines WAVE, 

RAYCON, REFRAC, CURVE, and HEIGHT are similar to those of their counterpart 
routines in Dobson's model (ref. 4). Subroutine DEPTH is an expanded version 
of its counterpart in Dobson's model, with additional logic for alternate 
topography approximation techniques and for determination of the required 
bathymetry data module. Subroutine WRITER stores computed parameters in arrays 
and prints these arrays after a complete ray path has been constructed. Sub- 
routine RAYINIT computes the spatial coordinates of the initial point on each 
ray in a particular set, based on the initial ray direction and a specified 
distance between adjacent rays and wave crests. Subroutine RAYOPT provides the 
user, if desired, with control over the final nearshore spacing between adjacent 
rays. This control option is also exercised if the user desires a plot of wave 
crests, which are constructed by segments in subroutine PLOTR. This subroutine 
also contains the logic for plotting both wave crests and wave -rays. Finally 
subroutine MATRIX performs matrix multiplications which are required in deter- 
mining coefficients in the bottom topography approximating polynomials. In 
addition, the program uses a number of system subroutines to satisfy various 
data handling and plotting requirements; these subroutines are listed in 
appendix A. 


Main Program WAVE 

The main program WAVE is used for input and output of parameters pertinent 
to the ray set being computed, including setting default values and checking 
for input errors. The random-access disk file named TAPE1 , which is used for 
input of bathymetry data modules, is initialized by using the system subroutine 
OPENMS. Initialization of parameters for each ray in the set is also performed 
in WAVE. If one of the plotting options is chosen, subroutine PLOTR is called 
to perform plot maintenance such as to open and close the plot file and to con- 
trol frame advances and positioning of frame origins. The flow diagram and 
listing of program WAVE follow. 
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PROGRAM WAVE ( INPUT^ OUTPUT»TAPE5«INPUT» TAPE 1 ) 


A PROGRAM TO CONSTRUCT REFRACTION DIAGRAMS AND COMPUTE WAVE 
HEIGHTS for WAVES MOVING INTO SHOALING WATER. 

♦♦♦♦ THIS PROGRAM REQUIRES A GRID OF BATHTMETRY DATA TO BE INPUT TO 
MEMORY. CREATE A RANDOM ACCESS DISK FILE (TAPED TO STORE THE 
DEPTH VALUES IN MODULE FORM. EACH MODULE IS TO CONTAIN A2 
SUB-ROWS AND 32 SUB-COLUMNS. SUB-COLUMNS 30>31* AND 32 OF A 
MODULE OVERLAY SUB-COLUMNS 1«2#AND 3 OF THE ADJACENT MODULE IN 
THE NEXT COLUMN. SUB-ROWS AO.Al.ANO A2 OF A MODULE OVERLAY 
SUB-ROWS 1»2/ AND 3 OF THE ADJACENT MODULE IN THE NEXT ROW. 
DEPTH AT POINTS IN THE ROWS AND COLUMNS WHICH LIE OUTSIDE THE 
REGIONAL MI BY MJ BATHYMETRY ARRAY ARE SET EQUAL TO THE DEPTHS 
ALONG THE BOUNDARIES OF THE ARRAY. THE MODULES ARE INDEXED BY 
ROWS BEGINNING IN THE COLUMN ADJACENT TO THE X-AXIS. 

♦*** NAMELIST INPUT PARAMETERS ♦♦♦♦* 

♦* BASIC DATA - NPUTl ♦* 

DCON - MULTIPLIER TO CONVERT DEPTH UNITS TO FEET FOR 

INTERNAL COMPUTATIONS 
(DEFAULT • 1.) 

DELTAS ■ MINIMUM STEP LENGTH ALONG RAY IN SHALLOW WATER IN 

GRID UNITS 

(DEFAULT • 0.01) 

GRID ■ NUMBER OF FEET PER GRID DIVISION 

(INPUT REQUIRED) 

GRINC > STEP LENGTH ALONG RAY IN DEEP WATER IN GRID UNITS 

(DEFAULT • 1.) 

lOUTUNT • OUTPUT UNITS - NAUTICAL MILES FOR SPATIAL 
COORDINATES X.Yj FDR ALL OTHER PARAMETERS# 

1 • U.S. CUSTOMARY UNITS 

2 • SI UNITS (DEFAULT) 

JPRTFRQ(l) • PRINT FREQUENCY FOR DEPTH/WL > 0.5 
(DEFAULT -20) 

JPRTFRQ(2) • print FREQUENCY FOR 0.25 < DEPTH/WL < 0.5 
(DEFAULT ■ 10) 

JPRTFR0(3) ■ PRINT FREQUENCY FOR DEPTH/WL < 0.25 
(DEFAULT • 5) 

KPLOT • PLOT OPTION 

0 ■ NO PLOT (DEFAULT) 

1 • PLOT WAVE RAYS WITH UNIFORM INITIAL SPACING 

2 " PLOT WAVE RAYS WITH CONTROLLED FINAL SPACING 

3 - PLOT WAVE CRESTS 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* * ** 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

« 

* 

* 

* 

* 


LIMNPT • MAXIrtUM NUMBER OF COMPUTATIONAL POINTS ALQN6 RAY 
(DEFAULT • 400) 

MI - TOTAL NUMBER OF ROWS IN REGIONAL BATHYMETRY ARRAY 

(INPUT REQUIRED) 

MJ • TOTAL NUMBER OF COLUMNS IN REGIONAL BATHYMETRY 

ARRAY 

(INPUT REQUIRED) 

NOMQD - NUMBER OF 42 BY 32 MODULES NEEDED TO STORE MI BY 

MJ ARRAY. MAXIMUM ■ 100 
(INPUT REQUIRED) 

NOROW - NUMBER OF MODULES ALONG X-AXIS IN NOHCD 

(INPUT REQUIRED) 

NOSETS - NUMBER OF SETS OF RAYS TO BE COMPUTED 
(DEFAULT • 1) 

PLTLNGX • LENGTH OF X-AXIS FOR GRAPHICS DISPLAY (INCHES) 
(DEFAULT - 10.0) 

PLTLNGY ■ LENGTH OF Y-AXIS FOR GRAPHICS DISPLAY (INCHES) 
(DEFAULT - 10.0) 

PLTSCAL • X-ANO Y-AXES SCALE FACTOR. CHANGE IN PHYSICAL X OR 
Y VALUE PER INCH OF PLOT AXES (NAUTICAL MILES PER 
INCH) 

(INPUT REQUIRED IF PLOTTING) 

TIDE ■ HEIGHT OF TIDE (FEET) 

(DEFAULT - 0.) 

SET DATA ■ NPUT2 ♦♦ 

♦ THE FOLLOWING SET OF DATA IS REQUIRED FOR ALL CASES. 

A . INITIAL DIRECTION OF RAY (DEGREES) 

(INPUT REQUIRED) 

HO ■ DEEP WATER WAVE HEIGHT (FEET) 

(INPUT REQUIRED) 

NCREST • SPACING OF CREST TIC-MARKS OR CREST CURVES IN 

INTEGER MULTIPLE OF THE PROGRAM COMPUTATIONAL TIME 
STEP 

(DEFAULT • 5) 

NSURFCE ■ BOTTOM TOPOGRAPHICAL APPROXIMATION TECHNIQUE 

1 - QUADRATIC LEAST SQUARES (DEFAULT) 

2 ■ CUBIC LEAST SQUARES 

3 • CONSTRAINED BICUBIC INTERPOLATION 

RAYSPC • INITIAL SPACING BETWEEN RAYS IN GRID UNITS 

(DEFAULT • 5.) 

T ■ WAVE PERIOD (SECONDS) 

(INPUT REQUIRED) 

* THE FOLLOWING PARAMETERS ARE REQUIRED ONLY IF KPLOT IS GREATER 

THAN OR EQUAL TO 2. 
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oistmax 


• MAXIHUM PERMISSIBLE FINAL SEPARATION DISTANCE * 

BETWEEN ADJACENT RAYS IN GRID UNITS * 

(DEFAULT • 999 .) ♦ 

DISTMIN • MINIMUM PERMISSIBLE FINAL SEPARATION DISTANCE * 

BETWEEN ADJACENT RAYS IN GRID UNITS * 

(DEFAULT ■ 25 .) ♦ 

KCREST - POINT BETWEEN ADJACENT RAYS AT WHICH TO BEGIN * 

PLOTTING THE CURVED CREST SEGMENT WHEN ANGLE ♦ 

BETWEEN adjoining FIRST-ORDER CREST SEGMENTS * 

IS LESS THAN RAYANG. MAXIMUM • 7 * 

(DEFAULT -1) * 

RAYANG ■ MINIMUM ACCEPTABLE ANGLE BETWEEN ADJOINING * 

FIRST-ORDER CREST SEGMENTS BELOW WHICH CREST * 

PLOTTING WILL BE MODIFIED (DEGREES) * 

(DEFAULT - 120.) ♦ 

RAYMAX ■ MAXIMUM PERMISSIBLE INITIAL RAY SPACING IN GRID * 

UNITS * 

(DEFAULT ■ 25 .) * 

RAYMIN - MINIMUM PERMISSIBLE INITIAL RAY SPACING IN GRID * 

UNITS * 

(DEFAULT ■ 0.01) * 

* 

*«*«*«* end of namelist input ******** * 

* 

** INPUT TITLE INFORMATION ♦♦ * 

* 

TITL ■ IDENTIFYING TITLE FOR EACH SET * 

(INPUT CARD FORMAT • 6A10) ♦ 

* 

**** DATA INPUT ORDER ♦♦♦♦* * 

* 

NPUTl/SET 1-TITL»NPUT2/S£T 2-TI TL . NPUT2/ETC * 

* 


*****t****************************************************************** 

* 

COMMON Bl» B2.C0.CXY.C PRIME. DCON» DELTAS# DEP(A2#32),DISTMAX/DI STM IN> 

♦ DRC.OTGR.DXY.GRIO. GRINC. HO. IBL K (1 01 ) » IGO. lOUTUNT. I RE F» ISTOP. JGO. 

♦ JPRTFRQ(3)»KCREST. KPLOT. L IMNPT. MI » MJ . NC REST. NFRST. NFRSTl. NL AST. 
♦NLASTl.NORAY.NaROW.NPT.NSURFCE.NTIC.NTICOFF.NWRITE.PDX.PDY.PD2X. 

♦ PD2 Y. PD2X Y. PL TL NGX. PLTLNG Y. PL TSCAL.R A YANG. RAYMAX. RAYMIN. RAYS PC. 
♦RCCO.RK.SK.Sl6.TI0E.TITL<8).V.WL»Wia»XO.XM.XPLOT{5OO)»XPLOTl(5OO)» 

♦ XPL0T2(500).XTEST. YO. YM. YPLOT ( 500 ). YPLOTl ( 500 ) . YPL0T2 ( 500 ). YTEST 

* 

N AME LI ST /NPUTl/DCON. DELTAS. GRID. GRINC. I OUTUNT.JPRTFRO.KPLOT. 
♦LIMNPT.MI.MJ.NOMOO.NOROW. NOSETS. PLTLN6 X.pl TLNGY. PL TS CAL. TIDE 
♦/NPUT2/ A. DISTMAX.O IS TM IN. HO.KC RES T.NCREST.NSURFCE. RAYANG. RAYMAX. 

♦ RAYMIN. RAYSPC. T 

* 

**** SET INPUT DEFAULT VALUES 
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* 


DCON ■ 1.0 
DELTAS - 0.01 
GRID • 0.0 
GRINC • 1.0 
lOUTUNT • 2 
IREPEAT • 0 
JPRTFRQ(l) • 20 
JPRTFR0(2I - 10 
JPRTFRQO) • 5 
KPLOT ■ 0 
LIMNPT ■ AOO 
MI ■ 0 
MJ - 0 
NOMOO - 0 
NOROW ■ 0 
NOSETS • 1 
PLTLN6X ■ 10. 

PLTLNGY ■ 10. 

PLTSCAL ■ 999. 

TIDE ■ 0. 

A ■ 999. 

DISTHAX • 999. 

OISTMIN ■ 25. 

HO ■ 999. 

KCREST • 1 
NCREST - 5 
NSURFCt • 1 
RAYANG ■ 120. 

RAYMAX - 25. 

RAYMIN • 0.01 
RAYSPC • 5. 

T • 999. 

* 

♦♦ READ BASIC DATA ♦♦ 

* 

READ (SfNPUTll 

* 

** CHECK FOR INPUT ERROR •• 

4 > 

IF (E0F(5)) 10,20 
10 PRINT 190 
GO TO 180 

20 IF (GRID .GT. O.OJ GO TO 30 
PRINT 200 
GO TO 180 

30 IF (MI .GT. 0 .AND. MJ .GT* 0 .AND* NOMDD .GT. 0 .AND, NOROW .6T« 
* 0 ) GO TO 90 
PRINT 210 
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<>0 IF (KPLOT .EO* 0 .OR. PLTSCAL «LT. 999.0) GO TO 50 
PRINT 220 
60 TO 180 

* 

** PRINT BASIC INPUT PARAMETERS ♦♦ 

♦ 

50 PRINT 230# DCON# OE LTAS# GRI 0#GRINC# lOUTUNT# J PRTFRQ# KPLOT# 

♦LIMNPT»MI#MJ»N0M00»NQR0W»N0SETS#PLTLNGX#PLTLNGY# PLTSCAL# TIDE 

A 

♦♦ OPEN MASS STORAGE DISK FILE (TAPED CONTAINING MODULAR BATHYMETRY DATA ♦♦ 

A 

NOMOD • NONOD+1 

CALL OPENMS ( 1# I8LK#N0M0D#0) 

A 

♦♦ OPEN PLOT FILE IF PLOT OPTION USED AA 

A 

IF (KPLOT .GT. 0) CALL PLQTR (1) 

A 

*+ GENERATE CURRENT RAY SET 

A 

DO 170 NOSET-l»NaS£TS 

A 

♦♦ READ title INFORMATION AA 

A 

READ 250» TITL 

A 

AA CHECK FOR TITLE CARD AA 
A 

IF (£0F(5)) 60# 70 
60 PRINT 260 
GO TO 180 
70 TITCl • 2H $ 

ENCODE (10#270#TITC2) TITL(l) 

IF (TITCl .NE. TITC2) GO TO 80 
PRINT 260 
GO TO 180 
80 CONTINUE 
A 

AA read input parameters for CURRENT SET ** 

A 

READ (5#NPUT2) 

A 

♦♦ CHECK FOR INPUT ERROR AA 

A 

IF (E0F(5)) 90#100 
90 PRINT 280 
GO TO 180 

100 IF (A .LT. 9999.0 .AND. HO .LT. 999.0 .AN[T. T .LT. 999.0) GO TO 110 
PRINT 290 


13 



6G TO 180 


*♦ INITIALIZE PARAMETERS FIXED FOR CURRENT SET ♦♦ 

* 

110 SIG ■ 6.28318531/T 
CC • 5.1204062+T 
WLO ■ CO*T 
DRC • 0.25*WL0 
OTGR ■ GRINC/CO 
AST - A 
DXY • 999.99 
ISTOP - 0 
NORAT ■ 1 
UNIT ■ OTGR*GRIO 

♦ 

♦♦ PLOT TITLE AND AXES INFORMATION IF PLOT OPTION USED ♦♦ 

* 

IF (KPLOT .GT. 0) CALL PLOTR (2) 

* 

♦♦ PRINT SPECIAL INPUT PARAMETERS IF PLOTTING WAVE CREST OR RAYS WITH ♦* 
♦♦ CONTROLLED FINAL SPACING ♦♦ 

IF (KPLOT .GE. 2) PRINT 2A0» DI STMAX# DI STMI N> KC RES T* RAYANG# 

♦ RAYMAX* RAYMIN 

♦ 

♦♦ FIND RAY STARTING POINT ♦♦ 

* 

120 CALL RAYINIT (X.Y.ASTJ 

IF (ISTOP .EQ. 1) GO TO 160 

* 

** INITIALIZE RAY PARAMETERS ** 

* 

A • AST 
CXY • CO 
NPT ■ 1 

NTICOFF ■ NCREST*(NTIC-1) 

NWRITE ■ 1 
WL - WLO 
B1«B2"RK»SK»1,0 
PDX-PDY-0. 

CALL DEPTH (X,YI 

* 

*♦ CALCULATE RAY PATH ♦♦ 

♦ 

CALL RAYCON (X.Y.A) 

* 

** SAVE INDEX OF FIRST AND LAST RAY POINT WITH RESPECT TO ♦♦ 

♦♦ REFERENCE CORNER ♦♦ 
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NFRST • NTICOFF+1 
NLAST ■ NTICOFF+NPT 

* 

♦♦ COMPUTE INITIAL RAY SPACING TO MEET FINAL CONSTRAINTS IF KPLOT *♦ 
♦♦ EQUALS 2 OR 3 *♦ 

* 

IF (KPLOT .GE. 2) CALL RAYOPT (IREPEAT) 

IF (IREPEAT .GT. 0) GO TO 120 

* 

** PRINT CONSTANT RAY PARAMETERS ** 

* 

PRINT 300# TITL#T»UNlT»RAYSPC#NCREST>NSURFCE»NORAY»NOSET 

* 

** PRINT OUTPUT HEADER INFORMATION ♦♦ 

* 

PRINT 310 

« 

♦♦ PRINT OUTPUT UNITS 

* 


IF (IOUTUNT-2) 130fl40 
130 PRINT 320 
GO TO 150 
140 PRINT 330 
150 CONTINUE 

♦ 

♦♦ PRINT COMPUTED PARAMETERS FOR CURRENT RAY ♦♦ 
CALL WRITER (X»Y*A#H) 

* 

♦♦ PLOT WAVE RAYS IF PLOT OPTION USED ♦♦ 


IF (KPLOT .GT. 0) CALL PLOTR (3> 
NORAY • NORAY+1 
GO TO 120 

160 PRINT 340# NOSET 
170 CONTINUE 

PRINT 350# NQSETS 
180 CONTINUE 

* 

** END PLOTTING IF PLOT OPTION USED ♦♦ 

« 

IF (KPLOT .GT. 0) CALL PLOTR (4) 
STOP 


190 FORMAT 
200 FORMAT 
210 FORMAT 
♦ ) 

220 FORMAT 


( 1H0.4INPUT 
(1H0»*INPUT 
( 1H0» *INPUT 

( 1H0»*INPUT 


ERROR 

■ 

NO 

ERROR 

m 

NO 

ERROR 

■ 

NO 

ERROR 

m 

NO 


DATA INPUT UNDER NPUTl*) 

DATA FOR GRI04) 

DATA FOR MI OR MJ OR NOMOD OR NOROW* 
DATA FOR PLTSCAL*) 
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230 FORMAT ( lHl>12X«*a ASIC PROGRAM PARAMETERS*// 

*2X>*DEPTH UNIT MULT.. 

*2X»*MIN STEP SIZE... 

*2X«*GRID SIZE 

*2X**0EEP WATER STEP SIZE 

*2X#*0UTPUT UNIT OPTION 

*2X«*PRINT FREO I (0/L>0.5) 

+2X**PRINT FREO 2 ( 0. 25 <0/L<0. 5 ) . . . 

♦2X#*PRINT FREO 3 <0/L<0.25> 

♦ 2X#*PL0T OPTION 

*2X«*MAX NO OF RAY POINTS 

*2X#*GRI0 LIMITS ABSCISSA 

*2X>* ORDINATE 

♦ 2X»*N0 OF DEPTH MODULES 

*2X#*N0 OF ROW MODULES IN NOMOD.... 

*2X#*N0 OF SETS OF RAYS 

*2Xf*PL0T X AXES LENGTH 

*2X.*PL0T Y AXES LENGTH.... 

*2X#*PL0T SCALE FACTOR 

*2X#*TI0E HEIGHT..... 

240 FORMAT <1H0»1X»*MAX FINAL RAY SPAC 

*2X>*MIN final ray SPACING 

♦2X#*INITIAL PT FOR CREST CURVE.... 

♦2X,*MIN ANGLE BETWEEN SEGMENTS.... 

*2X«*MAX INITIAL RAY SPACING 

■•■2X>*MIN INITIAL RAY SPACING 

250 FORMAT (8A10) 

260 FORMAT (1H0**INPUT ERROR - NO TITLE CARO*) 

270 FORMAT (A2) 

280 FORMAT (1H0«*INPUT ERROR - NO DATA INPUT UNDER 
290 format (1H0»*INPUT ERROR - NO DATA FOR A OR HO 


DC ON 

•*#F6.2/ 

DELTAS 

■*»F6.2/ 

GRID 

■*»F8.2/ 

GRINC 

•*»F6.2/ 

lOUTUNT 

■*»I2/ 

JPRTFRO(l) 

■*»I3/ 

JPRTFR0(2) 

-*.13/ 

JPRTFROO) 

•*.13/ 

KPLOT 

-*.12/ 

LIMNPT 

•*.14/ 

MI 

•*. 14/ 

MJ 

-*.14/ 

NOMOD 

•*.13/ 

NUROW 

•*.13/ 

NOSETS 

•*.13/ 

PLTLNGX 

•*.F6.2/ 

PLTLNGY 

•♦.F6.2/ 

PLTSCAL 

-*»F6.2/ 

TIDE 

•*.F6.2) 

ING 

..OISTMAX 

DISTMIN 

•*»F6.1/ 

KCREST 

•*. 14/ 

RAYANG 

•♦.F6.1/ 

RAYMAX 

•*.F6.1/ 

RAYMIN 

-*.F7.2) 


NPUT2*) 
OR T*) 


•♦.F6.1/ 


300 FORMAT (1H1,8A10/1X,*P£RIOD •♦.F6.2,* S EC . . .*»*TIME STEP -P.FT.Z. 
♦* SEC...*#*RAY SPACING -*» F5. 1» *». .CREST SPACING ■P.IB. 

**... BOTTOM APPROX CODE •*I2/1X.*RAY NO >*. I 3. *• • .SET NO ■*«I3) 


310 FORMAT (1H0»* POINT X 

*H SPEED REFRACTION SHOALING 

Y 

HEIGHT*) 

ANGLE 

DEPTH 

LEN6T 

320 FORMAT ( IIX. *( N.MI . ) 
♦) COEF COEF 

CN.MI.) 

(FT)*/) 

(DEG) 

(FT) 

(FT) 

(FT/SEC 

330 FORMAT ( 1 IX. * ( N. MI . ) 
♦ COEF COEF 

(N.MI.) 

(M)*/) 

(DEG) 

(M) 

(M) 

(M/SEC) 

340 FORMAT (1H0.*SET NO*. 
350 FORMAT (1H0#*ALL SETS 
END 

13.* ENDED. 
COMPLETED. 

RAYS LIE 
NUMBER OF 

OUTSIDE 
SETS -* 

BOUNDARY* 

.13) 

) 



Subroutine RAYINIT 


Subroutine RAYINIT computes the coordinates of the initial point on each 
ray in a particular set. It is assumed that the coastline and isobaths for the 
geographical region to be considered are aligned in a crude sense with the 
x-axis . Rays are to be initiated, depending on their initial direction, at a 
sequence of deepwater points along the seaward computational boundaries of the 
region (xQ,yj,,Xm), which are set arbitrarily a distance of 4 grid units inward 
from the limits of the bathymetry grid. (See fig. 1.) 

The initial ray in the set is found by forming a counterclockwise sequence 
of trial rays beginning at a reference corner of the region. For initial ray 
directions Oq > -90°, the reference corner lies at (xQ,y„]), whereas for 
ttg ^ -90°, it lies at (xnj,yjn). For olq = -180°, the initial ray is begun at 
the reference corner (Xm,ym), and no sequence of trial rays is formed. For 
all other initial directions, the counterclockwise sequence of trial rays 
(illustrated in fig. 1 for cxq > -90°) is formed, with each ray separated by 
a distance of RAYSPC grid units orthogonal to the initial ray direction. Steps 
are made inward along each trial ray in increments of DTIC = NCREST»GRINC 
until the test point along the ray lies within the appropriate seaward compu- 
tational boundary. If that point lies in deep water, the counterclockwise 
sequence is continued. If the first test point lying within the computational 
boundary is in intermediate water, the counterclockwise sequence is terminated 
and the previous ray in the sequence is chosen as the initial ray in the set. 
The sequence is also terminated if the entire trial ray lies outside the compu- 
tational boundaries or if the first intraboundary point along the trial ray 
lies on land. In these instances, the previous ray in the sequence is also 
chosen as the initial ray in the set. 

The coordinates of the initial point on each ray in the set are found by 
reversing the procedure and forming a clockwise sequence of rays beginning with 
the initial ray (illustrated in fig. 2 for Oq > -90°). The initial point 
along each ray is selected as that trial point advanced along the ray by an 
increment of DTIC and lying both within and closest to the appropriate seaward 
boundary. The clockwise sequence (and the ray set) is terminated when a trial 
ray is formed along which no appropriate initial point can be found. 

The flow diagram and listing of subroutine RAYINIT follow. 
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SUBROUTINE RAYINIT (X,Y#AST) 


t***t******* ************************** ******************************** 


* THIS SUBROUTINE COMPUTES THE COORDINATES OF THE INITIAL POINT ♦ 

* FOR EACH RAY IN A CLOCKWISE SEQUENCE ALONG THE SEAWARD * 

* BOUNDARIES OF THE REGION. EACH RAY IS SEPARATED BY RAYSPC GRID * 

* UNITS FROM THE PREVIOUS ONE. COMPUTATIONS BEGIN AT A POINT ♦ 

* ALONG THE RAY WHICH LIES IN DEEP WATER AND IS ADVANCED BY AN * 

* INTEGER MULTIPLE OF THE DEEP-WATER STEP LENGTH BETWEEN CREST ♦ 

* TIC MARKS(OTIC) FROM THE CREST WHICH WOULD PASS THROUGH THE ♦ 

* REFERENCE CORNER. ♦ 

^***t*t***************************************************************** 

COMMON Blf B2)C0 »CX Y>C PRIME .DC ON.OELTAS.de P( 4 2 >32) « 01 STMAX#DI STM INf 
>DRC. OTGR.DXY.GRID.gr INC. HO. IBLK (101). IGO. lOUTUNT. I REF. ISTOP. JGO. 

♦ J PRTF RQ ( 3 ). KCR E ST. KPLO T. L IMNPT. MI. MJ.NC RE ST. NFRST.NFRST1.nl AST. 
♦NLAST1.NORAY.NOROW.NPT.NSURFCE.NTIC.NTICOFF.NWRITE.PDX.PDY.P02X. 
>PD2Y.P02XY.PLTLNGX.PLTLNGY.PLTSCAL«RAYANG.RAYMAX.RAYM1N.RAYSPC. 
4^RCCO.RK.SK.SIG>TIOE.TITL(8).V. WL.WL0.X0.XM.XPL0T(500).XPL0T1 (300). 

♦ XPL0T2(500).XTEST, YO. YM. YPL OT ( 500 ). YPLOTl ( 500 ) . YPL0T2 ( 500 ). YTEST 

* 

IF (NORAY .EQ. 1) GO TO 10 
DXRAY— RAYSPC*SARAD 
DYRAY-RAYSPC*C ARAD 
GO TO no 

* INITIALIZATION SECTION 

* COMPUTE FIXED DEEP-WATER SEPARATION DISTANCES FOR RAYS AND CREST 

* TIC-MARKS. INITIALIZE COMPUTATIONAL BOUNDARIES. 

10 X0«4, 

Y0"4. 

XM«Ml-4 

YM-MJ-4 

ARAD-AST/57. 2957795 
SARAD-SIN( ARAD) 

CARAO-COS(ARAO) 

DXRAY»-RAYSPC*SARAD 

DYRAY-RAYSPCACARAO 

DTIC»NCREST*GRINC 

DXTIC-OTICACARAO 

DYTIC-DTIC*SARA0 

* 

* CHOOSE REFERENCE CORNER - (XO.YM) IF RAY ANGLE > -90 DEG. 

♦ (XM.YM) IF RAY ANGLE < -90 DEG. 
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IF (AST. LT. -89. 99) GO TO 20 

XC-XO 

YC-YM 

60 TO 30 

20 XC-XN 
YC-YH 

* LOCATE FIRST RAY IN SET 

* IF RAY ANGLE IS -180 DEG« FIRST RAY BEGINS AT REFERENCE CORNER. 

30 IF (AST. 6T. -179. 99) 60 TO ‘>0 
X-XTEST-XC 
Y-YTEST-YC 
NTIC«NC0RN»1 
RETURN 

* 

* ANGLE > -180 OEG. START SEARCH FOR INITIAL RAY AT REFERENCE CORNER. 

<»0 NTIC-1 
NCORN-0 
XTEST-XC 
YTEST-YC 

* 

* ADVANCE COUNTERCLOCKWISE ONE RAY AND HOVE TO CORRESPONDING TEST POINT. 

* 

50 XTEST-XTEST-DXRAY 
YTEST-YTEST-DYRAY 
NC0RN«NC0RN+1 

A 

♦ BASED ON AN6LE» TEST TO SEE IF CURRENT TEST POINT IS INSIDE REGION. 

♦ 

60 IF (AST. LT. -89. 99) 60 TO 70 
IF (XTEST.lt. XO) 60 TO 90 
IF (YTEST.LT.YO) GO TO ICO 
GO TO 80 

70 IF (YTEST.GT.YH+0.01) 60 TO 90 
IF (XTEST.LT.XO) GO TO 100 

* 

♦ INTERIOR TEST POINT FOUND. CHECK LOCAL DEPTH. 

* 

80 CALL DEPTH ( XTE ST. YTE ST ) 

♦ 

♦ IF DEEP WATER. ADVANCE COUNTERCLOCKWISE ONE RAY AND MOVE TO 

♦ CORRESPONDING TEST POINT. IF INTERMEDIATE WATER. THE PREVIOUS 

♦ RAY IS THE INITIAL RAY - GO TO STATEMENT 110 TO MOVE TO INITIAL RAY. 

♦ 

IF (DXY.LE.DRC) GO TO 110 
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GO TO 50 


* CURRENT TEST POINT OUTSIDE REGION. ADVANCE ONE TEST POINT 

* ALONG CURRENT RAY. 

* 

90 XTEST-XTEST+DXTIC 
YTEST-VTEST+DYTIC 
NTlC-MTIC+1 
GO TO 60 

♦ 

♦ ALL TEST POINTS ON TRIAL RAY OUTSIDE REGION, SET PREVIOUS RAY 

♦ TO THE INITIAL RAY. 

♦ 

100 CONTINUE 

* 

* ADVANCE CLOCKWISE ONE RAY AND MOVE TO CORRESPONDING TEST POINT* 

* 

110 XINIT-XTEST+DXRAY 
YINIT-YTEST+DYRAY 

♦ 

* IF THE RAY PASSING THROUGH THE REFERENCE CORNER HAS BEEN COMPUTED* 

* CRITERIA FOR ENDING RAY SET MUST BE TESTED. IF NOT* HAKE SURE THE 

* PRESENT TEST POINT IS THE FIRST INTERIOR TEST POINT. 

* 

IF (NORAY. GT. NCORN ) GO TO 160 

IF (XINIT.GT.XN^.OI.OR.YINIT.GT.YH^.OI) GO TO 160 
120 XTEMP"XINIT*DXT1C 
YTEMP-YINIT-OYTIC 
IF (AST, LT. -89, 99) GO TO 130 

IFIXTENP.LT.XO-O.Ol.OR.YTEMP.GT.YN+O.Ol) 60 TO 150 
GO TO lAO 

130 IF (XTEMP.GT.XM+.Ol.OR.YTEHP.GT.YM^.Ol) GO TO 150 

A 

♦ NEW INITIAL TEST POINT. RESET COORDINATES AND REDUCE VALUE OF TIC-NARK 

♦ COUNTER BY 1* GO BACK TO STATEMENT 120 AND TRY AGAIN, 

♦ 

lAO XINIT-XTENP 
YINIT-YTENP 
NTIC-NTIC-1 
GO TO 120 

* 

* CORRECT INITIAL TEST POINT FOUND. 

* 

150 X-XTEST-XINIT 
Y-YTEST-YINIT 
RETURN 

* 

♦ RAY PASSING THROUGH REFERENCE CORNER HAS ALREADY BEEN COMPUTED. 

* TEST VARIOUS CRITERIA FOR ENDING RAY SET, 
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160 IF (AST.LT.-.OI ) GO TO 170 

♦ 

♦ ZERO degree case. FINAL RAY ALREADY COMPUTED. 

* 

ISTQP-1 

RETURN 

170 IF (AST. LT. -90. 01) GO TO 210 

* 

* RAY ANGLE > -90 DEG. CHECK LOCATION OF PRESENT TEST POINT. 

* 

180 IF { YINIT.LE.YM) GO TO 190 

♦ 

* Y-COORDINATE OF CURRENT TEST POINT GREATER THAN YM. ADVANCE ONE 

* TEST POINT AND RECHECK LOCATION. 

A 

XINIT-XINIT+DXTIC 
YINIT-YINIT+DYTIC 
NTIC-NTIC+1 
GO TO 180 

190 IF (XINIT.lt. XM) GO TO 200 

* 

* X-COORDINATE OF CURRENT TEST POINT GREATER THAN XM. NO FURTHER RAYS 

* POSSIBLE. 

* 

ISTQP-1 

RETURN 

200 X-XTEST-XINIT 
Y-YTEST-YINIT 
RETURN 

♦ RAY ANGLE < -90 DEG. CHECK LOCATION OF CURRENT TEST POINT. 

* 

210 IF (XINIT.LE.XM+0.001) GO TO 220 

♦ 

* X-COOROINATE OF CURRENT TEST POINT GREATER THAN XM. ADVANCE ONE TEST POINT 

* AND RECHECK LOCATION. 

* 

XIN1T"XINIT*DXTIC 
YINIT-YINIT+DYTIC 
NTIC-NTIC+1 
GO TO 210 

220 IF(YINIT.LE.YO) GO TO 230 

* 

* INTERIOR TEST POINT FOUND. CHECK LOCAL DEPTH. 

* 

CALL DEPTH (XINIT. YINIT) 

IF (DXY.LE.DKC) GO TO 230 


♦ NEW RAY CAN BE COMPUTED 

♦ 

X-XTEST-XINIT 

Y-YTEST-YINIT 

RETURN 

♦ 

* INITIAL TEST POINT LIES EITHER OUTSIDE REGION OR IN INTERMEDIATE WATER 

* NO FURTHER RAYS POSSIBLE. 

* 

230 ISTOP-1 
RETURN 
END 
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Subroutine RAYCON 


Subroutine RAYCON controls the computation of each point along the wave 
ray. Points which lie in deep water are computed directly in RAYCON, while the 
routine calls subroutine REFRAC to compute ray points lying in shoaling water. 
After the ray point is computed, subroutine HEIGHT is called to compute the 
wave height and subroutine WRITER is called to store computed parameters for 
printing and plotting. Computation of the ray path is terminated if the ray 
goes outside the computational boundaries, if the shore is reached, or if the 
maximum number of ray computational points is exceeded. The program flag 
NWRITE designates the type of termination condition which is met. The differ- 
ent phases of the ray path are indicated by the program parameter IGO. When 
the ray is in deep water*, IGO = 1. When the ray reaches shoaling water, IGO 
is set to 2 in subroutine REFRAC. When a termination condition has been 
reached, IGO is set to 3 in subroutine WRITER. The flow diagram and listing 
of subroutine RAYCON follow. 
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SUBROUTINE RAYCON (X>Y>A) 


^t********************** *********************************************** 

* THIS SUBROUTINE ASSUMES DEEP WATER AND CALCULATES THE NEXT 

* POINT ON THE RAY. IF THE WAVE HAS ENTERED SHOALING WATER# THE 

* PREVIOUSLY CALCULATED POSITION OF THE NEXT POINT IS DISCARDED. 

* SUBROUTINE CURVE IS CALLED TO CALCULATE THE INITIAL RAY 

* CURVATURE# WHICH IS USED IN SUBROUTINE REFRAC TO CALCULATE THE 

* NEXT POINT ON THE RAY. 

^t********************** *************** ******************************** 

COMMON B1 # B2# CO#CX Y#CPRIME# DCON# DELTAS# DEP(A2# 32 )#DISTMAX#DISTMIN# 
♦DRC#DTGR#DXY*GRID»GRINC#HO#IBLK(101)#I60*IOUTUNT#IREF#ISTOP# J60# 
♦JPRTFRQ(31#KCREST#KPL0T#LIMNPT#MI»MJ»MCREST#NFRST#NFRST1#NLAST# 
♦NLASTl#NORAY#NOROW#NPT#NSURFCE#NTIC#NTlCOFF#NWRITE#PDX#PDY#PD2X# 
♦PD2Y#PD2XY#PLTLNGX#PLTLNGY,PLTSCAL#RAYANG#RAYMAX#RAYHIN#RAYSPC# 
4^RCCO#RK#SK>SIG«TIDE#TITL(8)#V#WL#WLO#XO#XM#XPLOT(500)#XPLOT1(500)# 
♦XPL0T2<500)#XTEST» YO#YM#YPLOT(500»#YPLOT1(500)#YPLOT2(500)#YTEST 

* 

♦♦ SET RAY INITIAL CONDITIONS ♦♦ 

* 

ANG ■ A/57.2957795 
COSA ■ COS(ANG) 

SINA > SIN(ANG) 

H - HO 
IGO • 1 

IF (NPT .EQ. 1) GQ TO 70 
10 PX • X 
PY ■ Y 

* 

** COMPUTE NEXT POINT ON RAY FOR DEEP WATER ♦* 

* 

DS • GRINC 
X « X*DS*COSA 
Y • Y+DS^SINA 

** TEST TO SEE IF RAY POINT OUTSIDE COMPUTATIONAL BOUNDARIES ♦♦ 

* 

IF (X .GE. XM .OR. X .LE. XO) GO TO 20 
IF (Y .LT. YM .AND. Y .GT. YO » GO TO 30 

* 

* OUTSIDE COMPUTATIONAL BOUNDARIES. SET FLAG AND EXIT ♦♦ 

* 

20 NWRITE ■ <t 
GO TO 70 


25 



** COMPUTE DEPTH ** 

* 

30 CALL DEPTH (X»Y) 

♦ 

♦♦ TEST TO SEE IF SHORE HAS BEEN REACHED ** 

* 

IF (OXY .6T. 0.) GO TO 40 

♦ 

♦♦ SHORE REACHED. SET FLAG AND EXIT *♦ 

* 

NWRITE • 3 
GO TO 70 

40 NPT - NPT ♦ 1 

* 

** TEST TO SEE IF MAX NO OF RAY POINTS REACHED ♦♦ 

* 

IF (NPT .LT. LIMNPT) GO TO 50 

* 

♦♦ MAX NO OF POINTS REACHED. SET FLAG ♦* 

♦ 

NWRITE - 5 

* 

** TEST TO SEE IF IN SHOALING WATER *♦ 

* 

50 IF (IGO .GE. 2) GO TO 60 
IF (OXY .GE. ORC) GO TO 70 

* 

*♦ IN SHOALING WATER. START REFRACTION CALCULATIONS *♦ 

* 

X • PX 
Y • PY 

♦ 

♦♦ COMPUTE INITIAL RAY CURVATURE ♦* 

* 

CALL CURVE (X»Y»ANG»FK) 

* 

♦♦ compute next POINT OS RAY FOR SHOALING WATER ** 

♦ 

60 CALL REFRAC (X/Y^ANG^FK) 

* 

♦* COMPUTE WAVE HEIGHT ♦♦ 

* 

CALL HEIGHT (ANG.H) 

* 

** SAVE PRINT AND PLOT DATA ♦♦ 

« 

70 A - 57.2957795*ANG 

CALL WRITER (X,Y,AfH) 

* 

** REPEAT FOR NEXT POINT ON RAY ** 

* 

IF (IGO-2) 10>40#80 
80 RETURN 
END 
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Subroutine REFRAC 


In subroutine REFRAC, equations (1) to (3) are solved iteratively to 
find the next point on the ray. If the curvature is oscillating between two 
values, the curvature is taken as the average of these two values. If the 
iteration process is converging very slowly or not at all, the ray is stopped. 
Subroutine CURVE is called to compute the curvature. Several ray termination 
conditions are checked in this subroutine. The ray is terminated if it goes 
outside the computational boundaries or if it reaches the shore. In addition, 
the ray is terminated if the wave speed has been reduced so that the incre- 
mental distance between successive ray points is less than the minimum step 
length allowed DELTAS which is an input parameter. The ray path phase flag IGO 
is set to 2 in this routine. The flow diagram and listing of subroutine REFRAC 
follow. 
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SUBROUTINE REFRAC (X>Y/A>FK) 


t*******************4******¥*************************^**************4* 


* THIS SUBROUTINE SOLVES THE REFRACTION EQUATIONS ITERATIVELY TO 

* FIND THE NEXT POINT ON THE RAY. IF THE ITERATIVE PROCESS IS 

* OSCILLATING BETWEEN TWO VALUES. THEY ARE AVERAGED. IF THE 

* PROCEDURE IS CONVERGING VERY SLOWLY, OR NOT AT ALL.THE RAY IS 

* STOPPED. 

t^L^^^itt**************************************************************** 

CONNON Bl. B2. CO. CXY. C PRIME. DCON. DELTAS. DEP( 4 2. 32 I.DISTNAX.DI STM IN, 
4 DRC, DTGR.DXY.gr ID. 6RINC.H0. IBLK( 101). 160. lOUTUNT. IREF. ISTOP. JGO. 
♦JPRTFR0(3),KCREST,KPL0T.LIMNPT,MI.MJ,NCREST,NFRST.NFRST1.NLAST. 
♦NLAST1.NORAY.NOROW.NPT.NSURFCE.NTIC.NTICOFF.NWRITE.PDX.PDY.P02X. 
♦P02Y.P02XY.PLTLNGX.PLTLNGY.PLTSCAL.RAYANG.RAYHAX.RAYMIN.RAYSPC. 
♦RCCO.RK.sk. SI G. TIDE.TITL( e) .V.WL.WLO.XO. XM. XPLOT (500 ).XPL0T1 ( 500) . 
♦ XPL0T2(500).XTEST. YO.YM. Y PLOT (500). YPLOTK 500). YPLCT2( 500). YTEST 

* 

NCUR ■ 1 

IF (IGO-2) 10.20.110 

♦♦ SET INITIAL CONDITIONS ♦♦ 

* 

10 FKM • FK 
IGO-2 

20 OS ■ CXYPDTGR 

* 

*♦ TEST TO SEE IF STEP SIZE TOO SMALL ♦♦ 

* 

IF (OS *GE« DELTAS) GO TO 30 

* 

♦♦ STEP SIZE TOO SMALL. ScT FLAG AND EXIT ** 

* 

NWRITE • 6 
GO TO 110 

A 

** COMPUTE TOLERANCE ON CONVERGENCE FOR CURVATURE ** 

* 

30 RESMAX > 0.00005/OS 

* 

** ITERATE FOR CURVATURE SOLUTION AND NEXT POINT ON RAY ♦* 

* 

40 DO 80 1-1.20 
DELA • FKM*OS 
AA - A^OELA 
AM - A^0.5A0ELA 
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XX ■ X+dS*CDS(AM) 
yy • y+DS*siN(AfO 


♦♦ TEST TO SEE IF RAT POINT OUTSIDE COMPUTATIONAL BOUNDARIES ♦* 

* 

IF (XX .CE. XM .OR. XX .LE. XO) GO TO 50 
IF (yy .LT. yM .and. yy .gt. yo) go to 60 

* 

♦* OUTSIDE COMPUTATIONAL BOUNDARIES. SET FLAG AND EXIT ♦♦ 

* 

50 NWRITE • 

GO TO 100 

♦ 

♦♦ COMPUTE RAT CURVATURE 

I* 

60 CALL CURVE ( XXj YY. AA» F KK ) 

* 

♦♦ TEST TO SEE IF SHORE HAS BEEN REACHED ** 

* 

IF (OXY .GT. 0.) GO TO 70 

* 

** SHORE REACHED. SET FLAG AND EXIT ** 

* 

NWRITE • 3 
GO TO 100 

70 IF (NCUR .EO. 2) GO TO 100 
FKM • 0.5*(FK*FKK) 

IF (I .EO. 1) GO TO 00 

* 

♦♦ TEST TO SEE IF SOLUTION HAS CONVERGED ♦♦ 

A 

IF (ABS(FKP-FKH) .LE. RESMAX) GO TO 100 
IF (I .EO. 18) FKie • FKM 
80 FKP ■ FKM 

* 

*♦ TEST TO SEE IF SOLUTION HAS CONVERGED ** 

* 

IF ( ABS(FKM-FK18) .LE. RESMAX) GO TO 90 

* 

** DIO NOT CONVERGE. SET FLAG AND EXIT *♦ 

* 

NWRITE • 2 
GO TO 100 

* 

** CURVATURE AVERAGE ♦* 

* 

90 FKM - 0.5*(FKM+FK18) 

NCUR • 2 
GO TO AO 

100 X - XX 

y ■ YY 

A ■ AA 
FK • FKK 
110 RETURN 
END 



Subroutine CURVE 


Subroutine CURVE initially calls subroutine DEPTH to determine whether 
the ray is in intermediate water or shallow water. It then sets the flag JGO 
to either 1 for intermediate water or 2 for shallow water. The local wave 
speed c and the rate of change of wave speed with depth c' are then com- 
puted according to the appropriate equations. By using these values, the 
ray curvature is computed. Several parameters required in subroutine HEIGHT 
for computation of the wave separation factor 3 are computed in this routine. 
The flow diagram and listing of subroutine CURVE follow. 
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SUBROUTINE CURVE (X»Y*A>FK) 

THIS SUBROUTINE CHECKS THE LOCAL DEPTH* CALCULATES THE WAVE 
SPEED FOR THE APPROPRIATE REGIME* COMPUTES THE RATE OF CHANGE 
OF SPEED WITH RESPECT TO DEPTH* AND CALCULATES THE CURVATURE 
AT THE SPECIFIED RAY POINT. 

^ 0 ***^)*********************************** ***************************** 

COMMON B1*B2*C0*CXY*CPRIME*DC0N*DELTAS*DEP(<>2*32)*DISTNAX*0ISTNIN* 
♦DRC*DTGR*DXY»GRID*GRINC*HO*IBLK(101)*IGO*IOUTUNT*IREF*ISTOP* JGO* 
4^JPRTFR0(3)*KCREST*KPL0T*LIMNPT*MI*MJ*NCREST*NFRST*NFRST1*NLAST* 

♦NL ASTI, NORAY* NOROW.NPT,MSURFCE.NTIC,NTICOFF*NWRITE*PDX*PDY*P02X* 
♦P02Y*PD2XY*PLTLNGX*PLTLNGY*PLTSCAL,RAYANG,RAYMAX*RAYMIN,RAYSPC* 
«RCCO,RK*SK*SIG*TIOE*TITL(8)*V*WL*WLO*XO*XM*XPLOT(500)*XPLOT1(500)* 
♦XPL0T2(500I,XTEST* YO*YM* YPLOT (5001 *YPLOTlt 500)* YPLOT2( 500)* YTEST 

IF (IGO .EQ. 1) 60 TO 20 

♦ COMPUTE DEPTH ** 

CALL DEPTH (X,Y) 

* TEST TO SEE IF IN INTERMEDIATE WATER ♦♦ 

IF (200.0*DXY .6T. WL ) 60 TO 20 
IN SHALLOW WATER ♦* 

• TEST TO SEE IF SHORE HAS BEEN REACHED *♦ 

IF (DXY . 6T. 0.) 60 TO 10 

♦ SHORE REACHED. EXIT SUBROUTINE *♦ 

GO TO 60 
10 J60 » 2 

* COMPUTE WAVE SPEED FOR SHALLOW WATER ♦♦ 

ARG ■ 32.1725*0XY 
CXY ■ SQRT(ARG) 

♦ COMPUTE CPRIME FOR SHALLOW WATER ♦* 

CPRIME - 16.08625/CXY 
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GO TO 50 


*♦ IK IKTERMEDIATE WATER ♦♦ 

♦ 

20 Cl ■ CXY 
JGO - 1 

* 

** ITERATE TO SOLVE FOR INTERMEDIATE WATER WAVE SPEED ♦♦ 

* 

DO 30 1-1*50 
ARG - (DXY*SIG)/C1 
CXY ■ C0*TANH(ARG) 

RESID • CXY - Cl 

IF (ABS(RESID) .LT. 0.00011 GO TO AO 
30 Cl - 0.5*(CXY+Cn 

* 

** COMPUTE INTERMEDIATE WATER CPRIME ♦* 

* 

40 Rcca • cxr/co 

SCMC - (l.-RCCO*RCCOJ*SIG 
V • SCMC4DXY+RCC04CXY 
CPRIME ■ CXY*SCMC/V 

* 

♦♦ COMPUTE RAY CURVATURE ♦* 

♦ 

50 FK - <SIN( A)*PDX-C0S(A)*PDYI*CPRIME/CXY 
60 RETURN 
END 
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Subroutine DEPTH 


Subroutine DEPTH first determines the bathymetry data module in which the 
current ray point resides. If the point lies outside the boundaries of the 
module currently residing in central memory, the required module is read into 
memory by using the system subroutine README. Subroutine DEPTH selects the 
subset of data points in the module that are to be used to compute the bottom 
topography approximation surface. It then, on option, computes coefficients 
to a least squares quadratic, a least squares cubic, or a constrained bicubic 
polynomial which is used to approximate the local topography. The depth at 
the current ray point is then computed by evaluating the approximating poly- 
nomial at that point. Partial derivatives required in subroutine HEIGHT for 
computing the wave separation factor 6 are also computed by evaluating 
analytical partial derivatives of the approximating polynomial. The input 
parameters TIDE and DCON are used to simulate a constant tide effect and to 
accommodate various units for the bathymetry data values. The flow diagram 
and listing of subroutine DEPTH follow. 
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SUBROUTINE DEPTH (X»Y) 


t**44******************************0********4***********4************* 


* THIS SUBROUTINE CALLS INTO CORE THE BATHYMETRY DATA MODULE IN * 

* WHICH THE CURRENT POINT RESIDES. IT THEN PITS A LEAST SQUARES * 

* QUADRATIC OR CUBIC POLYNOMIAL TO THE DEPTH VALUES AT 12 GRID * 

* POINTS SURROUNDING THE SPECIFIED POINT. ALTERNATIVELY. A * 

* CONSTRAINED BICUBIC POLYNOMIAL CAN BE USED TO FIT THE DATA BY ♦ 

* USING 16 GRID POINTS. THE DEPTH AND THE PARTIAL DERIVATIVES • 

* CAN THEN BE EVALUATED BY USING THE POLYNOMIAL APPROXIMATION. * 

4 4 

44444444444444 ** 4444444444444444444444444 *************** 4 *************** 

* 

COMMON B 1. B2. CO* CXY. CP RIME. DCON. DELTAS. DEP(A2>32)»DISTMAX#DI STM IN. 

♦ DRC. DTGR. DXY.GRID.gr INC. HO. IBLK( 101). IGO. lOUTUNT. IREF. ISTOP. JGO. 

J PRTF RQ( 3 ). KCRE ST. K PL OT.L I MNPT. Ml. MJ. NCREST.NFRST.NFRST1 .nl AST. 

♦NLAST1.NQRAY.NQR0W.NPT.NSURFCE.NTIC.NTIC0FF.NWR1TE.PDX.PDY.PD2X. 
♦PD2Y.PD2XY.PLTLNGX.PLTLNGY.PLTSCAL.RAYANG.RAYNAX.RAYMIN.RAYSPC. 
♦RCCO.RK.SK.SIG.TI0E.TITL(8).V.WL.WLO.XO.XM.XPL0T(50O).XPL0Tl(5OO). 
+XPL0T2(500).XTEST. YO.YM.YPLOTI 5001. YPLOTK 500). YPL0T2( 500). YTEST 

* 

DIMENSION C(16).0( 16).SX1 (6.12 ).SX2( 10.12).SX3(16.16) 

* 

* INPUT CONSTANT MATRIX FOR QUADRATIC LEAST SQUARES. SXl. BY COLUMNS . 

♦ 

DATA (SXl(I). I ■1.72)/. 306 6124 A. 2«. 05322967. -.125. .05263158. 

♦-.125. .2368 42 11.. 19677034. .105861 24. -.125. -.052631 5 8. -.125. 

♦.217703 3 5. 2 *.14413 876. -.12 5.. 05263158. -.12 5.. 2368 42 11 *.105861 24. 
♦.19677034. -.125. -.05263158. -.12 5. -.08492823*. 0903 11.. 0334928 2. 

♦.12 5.-. 15789474.0. .-.0514 3541. -.06758373*-. 03349282.. 125*. 15789474 
♦.0..-.05143541.-.0 3349282. -.06758 373*0... 15789474*. 125.-.08492823. 
♦.03349282.. 09031 1.0.. -.15789474*. 125.. 00598086.-. 18241627. 

♦. 12440191*. 125. -.15789474.0.. .13038278. -.3403 1101. -.12440191. .125. 
♦.15789474.0... 1303 8278. -.12440191. -.34031 101.0.. .15789474*. 125. 
♦.00598086*. 12440191. -.18241627*0.. -.15789474. .125/ 

♦ 

♦ INPUT CONSTANT MATRIX FOR CUBIC LEAST SQUARES. SX2. BY COLUMNS 

* 

DATA ( SX2( I ). I ■1.120)/. 3125.2*-. 625.-* 125*. 05263156. -.125. 4*. 25. 

♦.3125*. 625. -.625. -.125. -.052631 58. -.125. -.25*. 25. -.25*. 25*. 31 25. 
♦2*.625.-.125..05263158.-.125* 4*-.25. .3125.-.625. .625.-. 12 5. 

♦-.05263158. -.12 5.. 2 5. -.25.. 25. -.25. -.03125. -.0208 3 333*. 0625. .125. 
♦-.1578947*0.* .0833 3 33 3. -.25. 2*0.. -.031 2 5. -.02083 333.-. 062 5. .125. 

♦.157894 7.0... 08 333 3 33*. 25. 2*0.. -.031 25.-. 062 5. -.02083 333.0.. 

♦ .15 78947*. 125. 2*0... 2 5*. 08 3 333 33. -.0312 5*. 062 5. -.0208 3333.0.. 

♦-.1578947*. 125, 2*0 . 25, . 0833 333 3.-.03 125. . 020833 33.-.0625. 

♦. 12 5, -.1578947*0.. -.08 3333 33. .2 5. 2*0.. -.03 12 5*. 0208 3333*. 062 5, 

+. 125, .1578947*0. .-.08 3 33333.-. 25. 2*0. .-.03125*. 0625. .02083333. 

♦0... 1578947. .12 5. 2*0.,-. 2 5, -.08 33 333 3. -.03 12 5. -.062 5. .02083333. 

♦0.*-. 15 7894 7. .125. 2*0.. .2 5, -.08 33 333 3/ 

* 

* INPUT CONSTANT MATRIX FOR CONSTRAINED BICUBIC INTERPOLATION. SX3. BY COLUMNS 

♦ 

DATA (SX3( I).I-1.256)/1..2*0..-3..0..-3..2*.2*0..2..0..9,.-6..0.. 

♦— 6..4..0..1..0..-2..2*0..1..0.*— 3..0..2..6..— 4..0..— 3..2..2*0.*1.. 

♦ 2*0..— 2.,0., — 3..0. .1«.0.,6», — 3..2«. —4 •, 2., 4*0.. I*. 2*0*. 2 *— 2« . 0. . 

♦1..4..— 2..1.*— 2..1..3*0..3..2*0..— 2.,4*0..— 9..6..0..6.*- 4.,3*0.» 

♦— 1..2*0..1.,4*0..3..— 2..0.*- 3..2..7*0..3..3*0.*— 6..3.*— 2..4..— 2.. 

♦7*0..— 1..3*0..2*.-1..1..-2..1..11*0«.9..— 6..0..-6..4..11*0..— 3.. 
♦2..0..3.*- 2.*11*0. .— 3..3..0..2. .—2 ..11*0..1..— 1..0..— 1..1..5*0.. 

♦ 3..3*0..— 2..0..— 9. .6..0.* 6..— 4. .8*0. .3. .0..— 2..— 6.. 4..0.. 3.. — 2. . 

♦5*0..— 1..3*0.,1.,0.,3..— 3..0.*- 2..2..8*0.*— 1..0..1..2.*- 2..0.* 

♦-1..1./ 
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* 

IF (NSURFCE .EQ. 2) 60 TO 10 

* 

* PLACE LOCAL ORIGIN AT THE CENTRAL CELL CORNER CLOSEST TO X AND Y AXES FOR 

* CONSTRAINED BICUBIC INTERPOLATION AND QUADRATIC LEAST-SQUARES 

* 

XP-AMOD (X»l.) 

YP-AMOD (Y,l.) 

GO TO 20 

* 

* 

* PLACE LOCAL ORIGIN AT THE CENTER OF THE CENTRAL CELL FOR CUBIC LEAST-SQUARES 

♦ 


10 XP-AMOD (X,l.)-.5 
YP-ANOO (Y»l.)-.5 

* 

* CALCULATE REGIONAL INDICES OF NODAL POINT PI AND THE MODULE ROW AND 

* COLUMN NUMBERS USING GENERAL EQUATIONS. 

* 

20 Il-X+1. 

IM-Il/39+1 

Jl-Y+1. 

JM-Jl/29+1 

* 

IF (NPT.EO.l) GO TO 30 

* 

* IS CURRENT CELL SAME AS THE PREVIOUS ONE. IF S0> THERE IS NO NEED TO DO 

* ANYTHING AND THE PROGRAM SKIPS TO CALCULATION OF THE LOCAL OEPTH> DXY. 

* 

IF (I1.NE.HP) 60 TO 30 
IF (Jl.EQ.JlP) 60 TO 150 
30 IlP-Il 
JlP-Jl 


♦ CALCULATE LOCAL INDICES OF NODAL POINT PI USING GENERAL EQUATIONS. 

* 

ILl-MOO (Il>39) 

JL1«M00 (J1.29) 

IF (ILl.GT.l) GO TO 40 

* 

* IF ILl IS LESS THAN 2» HE MUST USE THE ALTERNATE EQUATIONS FOR MODULE 

* ROW NUMBER AND LOCAL ROW INDEX. THUS# SET IMADD>-1 AND ILlADD-39. 

* 

IMADD—1 
ILlADD-39 
GO TO 50 
40 IMAOO-0 
ILlADD-0 

50 IF (JLl.GT.l) GO TO 60 

* 

* IF JLl IS LESS THAN Z, WE MUST USE THE ALTERNATE EQUATIONS FOR MODULE 

* COLUMN NUMBER AND LOCAL COLUMN INDEX. THUS# SET JMAOO— 1 AND JLlADD-29. 

* 

JMADD — 1 
JLlADD-29 
GO TO 70 
60 JMAOD-0 
JLlADD-0 
70 IM-IM+IMADO 
JM-JM-I-JMADO 
I"IL1*IL1ADD 
J-JL14-JL1A0D 
INDEX-! JM-1)*N0R0W*IM 
IF (NPT.EO.l) GO TO 80 
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♦ IF INDEX • INOEXP, THE CURRENT NODULE IS THE SAHE AS THE PREVIOUS ONE 

♦ AND THERE IS NO NEED TO READ THE R ANDOM-ACCESS FILE. SKIP TO FILLING 

* THE D ARRAY. 

* 

IF (INDEX. EQ. INOEXP) GO TO 90 
80 INDEXP-INOEX 

* 

* READ MASS STORAGE DISK FILE (TAPED CONTAINING BATHYMETRY DATA MODULE 

* 

CALL READMS ( If DEP » INDEX ) 

* 

90 IF (NSURFCE.E0.3) GO TO 100 

* 

* FILL VECTOR FOR QUADRATIC OR CUBIC LEAST SQUARES APPROXIMATION 

* 

D(l) ■ DEP(DJ) 

D(2) • DEP(Dl.J) 

D(3) ■ DEP(I*1»J+1) 

D(9) > DEPdiJ^l) 

0(5) ■ DEP(I+2,J) 

D(6) - DEP(D2fJ4^1) 

D(7) • DEP(Dl»J+2) 

D(8) - 0EP(IfJ>2) 

D(9) - 0EP(I-1»J*1) 

D(10) ■ OEP(I-DJ) 

D(1D • OEP(I.J-l) 

D<12) ■ DEP(D1»J-1) 

GO TO 110 

* 

* FILL VECTOR FOR CONSTRAINED BICUBIC INTERPOLATION 

* 

100 D(l)*DEP(If J) 

0(2)-0.5*<OEP(I+1» J)-0£P< I-l# J)l 
0(3)-0.5*(DEP(l> J+D-0EP(D J-ID 

0(9)«0.25*(DEP(I*1» J+1)-DEP(I+1» J-l)-0EP(I-lf J+l)+0EP(I-lf J-1) ) 
D(5)-0EP(I*1»J ) 

D(6)>0.5*(DEP(D2f J)-0EP(If J) ) 

D(7)-0.5*(DEP(D1. J*1)-DEP( Dlf J-l)l 

0(8)-0.25*(OEP(I+2» J + D-DEP{I + 2»J-l)-DEP(If J + l)+DEP(If J-D) 
D(9)«DEP(I+1. J*l) 

0(10)»0.5*(DEP( 1+2* J + 1)-DEP(I» J-H) ) 

D(11)»0.5*(0EP( I + l» J*2)-DEP(D1» J)) 

0 (12)»0.2 5*(DEP(D2, J+2)-DEP( D2» J)-OEP(If J + 2)*DEP(I» J) ) 
D(13)>DEP(If J«'l) 

D(1<.)-0.5*(0EP(I41» J*1)-0EP(I-1» J+ll ) 
0(15)»0.5*(0EP(I.J+2)-DEP(I»J) ) 

D(16)«0.25*(DEP(I+1» J+2)-DEP(I*l» J)-0EP(I-1» J*2)*DEP( I-lf J) ) 

♦ 

110 IF (NSURFCE-2) 120fl30#190 

* 

* FIND COEFFICIENTS OF QUADRATIC LEAST SQUARES POLYNOMIAL BY 

* MATRIX MULTIPLICATION. CC] ■ CSXIHOI 

* 

120 CALL MATRIX ( 6# 12> SXD Of C ) 

60 TO 150 

* 

* FIND COEFFICIENTS OF CUBIC LEAST SQUARES POLYNOMIAL BY 

♦ MATRIX MULTIPLICATION. CC] • CSX21CD1 



130 CALL MATRIX ( 10» 12» SX2* D>C ) 

GO TO 150 

* 

♦ FIND COEFFICIENTS OF CONSTRAINED BICUBIC INTERPOLATING POLYNOMIAL BY 

♦ MATRIX MULTIPLICATION. CCl • CSX3]tOJ 

♦ 

1<*0 CALL MATRIX ( 16» 16» SX3j D* C ) 


* EVALUATE DEPTH BY USING POLYNOMIAL APPROXIMATION 

* 

150 DXY-C(1)+(C(3)+(C( 6)*C(10J*YP)*YP)*YP+(C(2)+(C(5)*(C(9)*C(11)AYP)* 
♦YP)*YP)*XP+(C(4)+(C(8I+(C(12>*C{13»*YP)*YP)*YP)*XP**2+(C(7)+ 
•KC(1‘») + (C(15)+C(16)*YP)*YP)*YP)*XP**3 

A 

♦ COMPUTE PARTIAL DERIVATIVES OF LOCAL DEPTH WITH RESPECT TO SPATIAL 

♦ COORDINATES FROM POLYNOMIAL APPROXIMATION 

A 

P0X-C(2)a(C(5)a(C(9)+C(11»ayp|AYPIAYP+2,a(C(A)a(C(8)a(C(12)aC(13)* 

aYP)ayP)ayp)axPa3.axpaa2A(C(7)a(C(1«)a(C(15)aC(16)AYP)AYP)AYP) 

PDY-C (3)a(C(5)a<C(8»aC(1<.)*XP>AXP)AXP+2.AYP*(C(6)a(C(9)+(C(12)a 
♦C(15)*XP)*XP)AXP)+3.*YPAA2A(C(10)+(C(11)+(C(13)aC( 16)AXP)AXP)AXP» 
PD2X«2.*(C(4)+(C(8)A(C(l2)ACtl3 )AYP) AYP)AYP)a6.AXPACC(7)4(C(1A)+ 
♦<C(15)*C(16)ayP)AYP)ayP» 

P02Y-2.A(C(6)+(C(9 )+(C (12 )*C (15 ) AXP ) AXP ) AXP ) aB.* YP A (C ( 10 ) A(C ( 11 ) ♦ 
♦(C(13)AC(16)AXP)AXP)AXP) 

PD2XY«C(5)a2.*(C(9)*YP+C{8)*XP)a3.*(C(11)*YPAA2+C(1A)AXPAA2)a 
♦4.AC(12)*XPAYP+6.AXPAYPA(C(13J AYPaC(15)axP»*9.AC(16)AXPAA2AYPAA2 

A 

♦ CHANGE DEPTH UNITS BY DCON 

A 

DXY-OXYADCON+TIDE 

PDX-PDXAOCON 

PDY-PDYAOCON 

PD2X-PD2XMDC0N 

PD2Y-PD2YADC0N 

P02XY-P02XYADC0N 

A 

RETURN 

END 
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Subroutine HEIGHT 


Subroutine HEIGHT calculates the local shoaling coefficient (eq. (7)) and 
then computes the refraction coefficient (eq. (8)) by using the value of 8 
calculated at the previous ray point. These coefficients are then used to 
compute the local wave height. Fox's finite-difference solution for wave 
intensity is solved to give a new value of 8. The flow diagram and listing 
of subroutine HEIGHT follow. 
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SUBROUTINE HEIGHT (A>H) 


^4tt******************************************************************* 

THIS SUBROUTINE CALCULATES THE SHOALING COEFFICIENT^ SK« AND * 
COMPUTES THE REFRACTION COEFFICIENT^ RK, USING THE VALUE OF ♦ 

BETA AT THE PREVIOUS POINT ON THE RAY. THESE ARE THEN USED TO * 
CALCULATE THE LOCAL WAVE HEIGHT. THEN THE FINITE-DIFFERENCE ♦ 

FORM OF THE EQUATION OF WAVE INTENSITY IS SOLVED TO GIVE A NEW * 
VALUE OF BETA. * 

^tt******************************************************************** 

COMMON B1.B2jC 0.CXY.CPRIME.DC0N.DELTAS.DEP(<»2.32).DISTHAX.DISTMIN. 
•^ORC>OTGR#DXY#GRID*GRINC«HO>I8LK<101)>IGQ#IOUTUNT#IREF#1STOP# JGO# 
•^JPRTFRQ(3)«KCREST.KPL0T«LIMNPT»MI*MJ*NCREST>NFRST#NFRST1#NLAST# 
>NLASTl.NORAY.NORaw.NPT>NSURFCE>NTIC.NTICOFF>NWRITE#POX#PDY*PD2X# 
tPD2Y#PD2XY»PLTLN6X>PLTLN6Y.PLTSCAL,RAYANG.RAYMAX.RAYMIN»RAYSPC» 
•^RCCO.RK.SK>SIG#TI0E#TITL(8)«V«WL>WLQ>XO>XM*XPLOT(50O)#XPL0Tl(5OO)> 
♦ XPL0T2(500)»XTEST» YO. YM, YPLOTI 500 ) » Y PLOTl ( 500 ) # YPL0T2 ( 500 )# YTEST 

IF (160 «LT. 2) GO TO 30 
WL ■ WLO^RCCO 
GN - 12.5663706*0XY/WL 
CG ■ (1.0*GN/SINH( 6NI)*CXY 
IF (CG .LT. 0.) GO TO 30 
SK • SQRT(CO/CG) 

♦ COMPUTE REFRACTION COEF USING MUMK-ARTHUR EON ♦♦ 

RK ■ ABS(1./B2) 

RK ■ SQRT(RK) 

'• COMPUTE WAVE HEIGHT ♦♦ 

r 

H • HO*SK*RK 

I 

TEST TO SEE IF IN SHALLOW OR INTERMEDIATE WATER ** 

I 

IF (JGO .EQ. 2) GO TO 10 

[ 

'♦ INTERMEDIATE WATER VALUE FOR U ♦♦ 

( 

U • -2.*SIG*RCC0*CXY/(V*V) 

GO TO 20 

SHALLOW WATER VALUE FOR U ** 

10 U • -0.5/0XY 

'* COMPUTE RAY SEPARATION FACTOR ♦♦ 

20 COSA > COS(A) 

SINA • SIN(A) 

P ■ -(C0SA*PDX+SINA*P0Y)*CPRIME*DTGR*2. 

0 • ( (PD2X*U*PDX*P0X)*SINA*SINA-(PD2XY+U*PDX*PDY)*2.*SINA*C0SA+ 
♦(P02Y+U*PDY*PDY)*C0SA*C0SA)*CPRIME*CXY*0T6R*DTGR*2. 

B3 • ((P-2.)«B1«(4.-0)*B21/(P4^2.) 

B1 > B2 
B2 - B3 
30 RETURN 
END 



Subroutine WRITER 


Subroutine WRITER is used for storing computed parameters in appropriate 
arrays and then for printing the same parameters after the complete ray path 
has been constructed. The output units for spatial coordinates (x,y) are in 
nautical miles. The output units for all other computed parameters can be in 
either U.S. Customary Units or the International System of Units (SI), as 
specified by the parameter lOUTUNT. The first and last points on a ray are 
always printed, and intermediate points are printed according to the print 
frequencies specified by the input array JPRTFRQ. Position data are saved in 
this routine to be used in the plotting subroutine PLOTR. The flag IGO is set 
to 3 when the ray has reached a termination condition, and a message specify-' 
ing the particular termination condition is printed. The flow diagram and 
listing of subroutine WRITER follow. 
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SUBROUTINE WRITER (X«Y»A#H) 


*^^t******************** *********************************************** 

* 

* THIS SUBROUTINE STORES AND PRINTS COMPUTED RAY PARAMETERS AND 

* THE RAY TERMINATION CONDITION. X AND Y COORDINATES ARE SAVED 

* IN THIS SUBROUTINE FOR EVENTUAL PLOTTING. 

* 

4^^^4t^^t************************************************************** 

* 

COMMON Bl>B2>COfCXYfCPRIME*DCON»DELTAS»DEP(^2«32)#DISTMAX»DISTMIN« 
4^DRC«DTGR«DXY»GRIDf GRINC»H0«IBLK(101) .IGO. lOUTUNT.IREF.ISTOP. JGO* 
♦JPRTFRQ(3)»KCREST»KPL0T#LINNPT>MI#MJ#NCREST>NFRST/NFRST1»NLAST» 
4^NLASTl«NORAYfNORaw>NPT>NSURFCE>NTIC>NTICOFF>NWRITE.PDX.POY.PD2X. 
«PD2Y*PD2XY*PLTLNGX«PLTLNGY>PLTSCAL>RAYANG#RAYMAX»RAYMIN#RAYSPC# 
4^RCCO«RK>SK»SIG>TI0E«TITL(a)>V#WL»WLO«XO*XM>XPLQT(&00)#XPL0Tl(5OO)# 
+ XPLQT2(500).XT£ST» YO, YM, YPLOT ( 500 ) » Y PLOTl ( 500 ) » YPL0T2 ( 500 )» YTEST 

* 

DIMENSION NPT1(400)>X1(AOO).Y1(<>00). Al(<>OO)>DXYl('<r0O>tWLl('«OO). 

♦ CXY1(400)>RK1(<*00) »SK1(400)»H1(<.00) 

* 

IF (IGQ .EQ. 3) GO TO 30 
IF (NWRITE .GE. 2) IG0>3 

* 

** SET OUTPUT UNITS *♦ 

* 

IF (NPT .GT. 1) GO TO 10 
NC ■ 0 
AMULT • 1, 

XYMULT • 1./6076.1155 

IF (lOUTUNT .EQ. 2) AMULT • 0.30A8 

* 

*♦ SET PRINT FREQUENCY ♦♦ 

* 

10 IPRINT • JPRTFRQ(3) 

IF (OXY/WL .GT. 0.25 .AND. DXY/WL .LT. 0.5) IPRINT > JPRTFR0(2) 

IF (DXY/WL .GE. 0.5) IPRINT • JPRTFRO(l) 

♦ 

♦♦ SAVE X AND Y RAY COORDINATES FOR PLOTTING WHICH ARE INDEX ♦♦ 

*♦ WITH RESPECT TO REFERENCE CORNER *♦ 

* 

XPLOT(NPT*NTICOFF) ■ X 
YPLOT(NPT+NTICOFF) • Y 

* 

** SAVE FIRST AND LAST POINT ♦♦ 

* 

IF (NPT .EQ. 1 .OR. NWRITE .Gc. 2) GO TO 20 



♦♦ PRINT OUTPUT CONTROL ♦* 

* 

IF (NOD (NPTfIPRINT) .NE. 0) GO TO 50 

* 

** SAVE COMPUTED RAY PARAMETERS ** 

* 

20 NC • NC ♦! 

NPTKNC) • NPT 
XKNC) • XYMULT*GRID*X 
Yl(NC) ■ XYMULT*GRIO*Y 
AKNC) - A 

DXYl(NC) ■ AMULT*OXY 
WLl(NC) * AMULTPUL 
CXYKNC) • AMULT*CXY 
RKl(NC) > RK 
SKKNC) - SK 
HKNC) ■ ANULT+H 
GO TO 50 

« 

♦♦ PRINT STORED ARRAY OF COMPUTED RAY PARAMETERS 

* 

30 DO AO I-1«NC 

AO PRINT 60 «NPTl(I)*Xl(I)«Yl(n«Al(n#DXYl(I)#WLl(I)/CXYl(I), 
^RKl(I)>SKl(I)>Hl(n 

* 

*♦ PRINT RAY TERMINATION CONDITIONS ♦♦ 

* 

IF (NWRITE .EQ. 2) PRINT 70 
IF (NWRITE «EQ* 3) PRINT 80 
IF (NWRITE ,EQ. A) PRINT 90 
IF (NWRITE *E0. 5) PRINT 100 
IF (NWRITE «EQ« 6) PRINT 110 
50 RETURN 
A 

60 FORMAT (3X#IA#2F10.2f2X/2F9.2#2F8.2#F10.3#F11.3#F9«2) 

70 FORMAT (IHO/A raY STOPPED* NO CONVERGENCE FOR CURVATURE*) 

80 FORMAT (1H0»* RAY STOPPED# REACHED SHORE*) 

90 FORMAT (IHO#* RAY STOPPED# REACHED BOUNDARY*) 

100 format (1H0»* ray STOPPED# NUMBER OF POINTS EXCEEDS MAXIMUM*) 

110 FORMAT (IHO#* RAY STOPPED# INCREMENTAL DISTANCE ALONG RAY LESS TH 
♦AN MINIMUM*) 

END 
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Subroutine PLOTR 


Subroutine PLOTR is used for segment-by-segment construction of wave 
crests and for plotting the ray diagrams. The logic for constructing the 
crest patterns is detailed in appendix B. The plotting mode is controlled 
by the input parameter KPLOT. Axes, labels, and notations are the same for 
both the ray and crest modes of plotting. The input parameters PLTLNGX and 
PLTLNGY control the size of the plotting area, and the input parameter PLTSCAL 
scales the data to fit this area. Actual plotting is performed by using system 
subroutines which are described in appendix A. The flow diagram and listing 
of subroutine PLOTR follow. 
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SUBROUTINE PLOTR (KOPT) 


* * 

* THIS SUBROUTINE PLOTS RAY AND CREST DATA, ♦ 

* * 
^,^^*^*^*t*************************************************************** 

COMHON Bl,B2/CO,CXy,CPRINE,OCON,DELrAS,DEP(42#32)#OISTNAX#OISTNIN> 
-^DRC*0TGR,DXY«GRI0,GRINC,H0#IBLK(101) , IG0> lOUTUNT, IREF, ISTOP, JGO, 
4^JPRTFRQ(3),KCREST,KPL0T«LINNPT>HI*NJ>NCREST«NFRST*NFRST1»NLAST» 
•^NLAST1*N0RAY,N0R0W*NPT,NSURFCE>NTIC*NTIC0FF>NURITE#PDX#PDY#P02X# 
♦PD2Y>PD2XY>PLTLNGX,PLTLNGY,PLTSCAL»RAYANG,RAYHAX,RAYNIN,RAYSPC» 
4-RCCQ,RK,SK>SIG,TI0E,TITL(8>«V,WL, WLQ« XO*XM« X PLOT ( 500) t XPLOTl ( 500 ) , 
♦XPLQT2(500),XT£ST. YO, YN,YPLOT(500),YPLOT1(500)#YPLOT2(500)#YTEST 

* 

DIMENSION XPL(<*00),YPL(400) 

* 

GO TO (10, 30,100, 230), KOPT 

* 

♦♦ OPEN PLOT FILE AND SET PLOT SPEED *♦ 

* 

10 CALL PSEUDO 
CALL LEROY 
EPS - 0.001 
DO 20 1-1,500 

20 XPLQTd )-YPLOT(l)-XPLOTl(I)-YPLOTltI )-XPL0T2( I)-YPL0T2(I)-0,0 
GO TO 240 

♦ 

*♦ ADVANCE PLOT FRAME AND SET CONSTANTS ♦♦ 

♦ 

30 CALL NFRAME 

AN-GRlO/6076.1155 
COSRAY»COS(RAYANG/ 57, 2957795) 

* 

** SET ORIGIN *♦ 

4 > 

CALL CALPLT (3. ,1.5, >3) 

* 

** PLOT TITLE INFORMATION ♦♦ 

* 

CALL NOTATE (-1,5, 1.0, 0,15, TITL, 90., 60) 

IF (NSURFCE-2) 40,50,60 

40 CALL NOTATE (-1 .2, 3. 0, 0.1 25,2 3H0UADR ATIC LEAST SQUARES, 90. ,23) 

GO TO 70 

50 CALL NOTATE (-1.2, 3.0,0. 125,19HCUBIC LEAST SQUARES, 90., 19) 

60 TO 70 

60 CALL NOTATE (-1.2, 3.0,0. 125,33HC0NSTRAINED BICUBIC INTERPOLATION, 
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4-90. » 33) 


♦♦ PLOT AND LABEL AXES *♦ 

70 Yl—0.069 

TMAJ-IO./PLTSCAL 

LMT«PLTLNGY*PLTSCAL4lO. 

00 80 I"10»LnT,10 
FPN-I-10 

CALL NUMBEtl (-0. 15 , Y 1* 0. 1 2. FPN» 90 .,-1 ) 

80 Yl-Yl+TMAJ 
Xl-0.06 

LHT-PLTLNGX+PLTSCAL+10. 

DO 90 I»10»LMT»10 
FPN-I-10 

CALL NUMBER ( Xl»-0 . A, 0 .12» FPN» 90. »-l ) 

90 Xl-Xl+TMAJ 

YH-PLTLNGY/2.0-0.26 

CALL NOTATE (-0 .6. YH. 0.1 5 » AHY. NM, 90.# A) 

XH-PLTLNGX/2.040.08 

CALL NOTATE {XM»-1. 20# 0. 1 5# AHX# NH# 90. # A ) 

CALL AXES (0.#0.«0.#PLTLNGX#0.#PLTSCAL#TNAJ#PLTSCAL#1H #0.#1) 
CALL AXES (0.#0.#90.#PLTLNGY#0.#PLTSCAL#TNAJ#PLTSCAL#1H #0.#-1) 
CALL CALPLT (0.,PLTLNGY»3 ) 

CALL CALPLT ( PLTLNGX, PLTLNGY# 2 ) 

CALL CALPLT (PLTLNGX#0.# 2) 

GO TO 2A0 

100 IF (KPLOT .£0. 3) GO TO lAO 

A 

*♦ PLOT RAY ♦♦ 

A 

NS-0 

DO 110 I-NFRST*NLAST 
NS-NS+1 

XPL(NS)«AN*XPLOT(I )/PLTSCAL 
110 YPL(NS)«AMPYPLOT(I)/PLTSCAL 

IF (MOO (NORAY.S) .NE. 0) GO TO 120 

A 

AA PUT RAY NUMBER ON PLOT AA 

A 

XNO - XPL(1)A0.03 
YNO ■ YPL(l)+0.03 

CALL NUMBER ( XNO, YNO# 0,07. FLOAT (MORAY )# 90. #-l ) 

A 

AA PUT A SYMBOL AT EVERY NCREST POINT ON RAY AA 

A 

120 DO 130 I-1#NS#NCREST 

130 CALL PNTPLT ( XP L ( I ) # YPL ( 1 1 #-22# 1 ) 

CALL DRAW (XPL#YPL#NS) 
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60 TO 190 

140 IF (NORAY .LT. 3) GO TO 190 

♦♦ OETERHINE FIRST ANO LAST RAY POINT WITH RESPECT TO REFERENCE ♦♦ 
♦♦ CORNER FOR PLOTTING CREST DATA ♦♦ 

* 

HSl ■ NFRST 

IF (NFRSTl .GT. NFRST) HSl-NFRSTl 

IF (NFRST2 .GT. NFRST .AND. NFRST2 .GT. NFRSTl) HS1-NFRST2 
HS2 • NLAST 

IF (NLASTl .LT. NLAST) MS2-NLAST1 

MSl-MSl-H 

MS2-MS2-1 

* 

** PLOT CREST DATA ♦♦ 

DO 180 I>NS1«MS2.NCREST 

* 

** CALCULATE VECTOR COMPONENTS OF RAY SLOPES AND FIRST-ORDER ♦♦ 

*♦ CREST SEGMENTS *♦ 

* 

DXS1»XPL0T1(I+1)-XPLQT1(I-1) 

D YS !■ YP LOTI ( 1*1 )-YPLOTl( 1-1) 

IF (A9S(DYS1) .LT. EPS) 0YS1»SI6N(EPS»0YS1 ) 


SL1--DXS1/DYS1 
0XPT-XPL0T(I+1 )-XPLOT(l-l) 

DYPT«YPLOT( I+1)-YPL0T( I-l) 

IF (ABS(DYPT) .LT. EPS) DYPT-S I 6N ( EP S» D YPT ) 

SL2—DXPT/DYPT 

OXPS-XPLOT(I)-XPLOTKI) 

IF (ABS(DXPS) .LT. EPS) DXPS-S IGN ( EPS. DXPS ) 

DYPS-YPLOTd )-YPLOTl(I ) 

D1«SC)RT(DXPS**2 + 0YPSPA2) 

IF (01 .LT. EPS) 01-EPS 
IF (I .GT. NLAST2) GO TO 160 
DXS12-XPL0T2(I )-XPLOTl(I ) 

0YS12«YPL0T2(I )-YPLOTl (I ) 

D2«SQRT(0XS12**2+0YS12**2) 

IF (D2 .LT. EPS) 02-EPS 

* 

CALCULATE ANGLE BETWEEN ADJOINING FIRST-ORDER CREST SEGMENTS ♦* 

THETA-(0XPS*0XS12t0YPSADYS12)/ (D1*D2) 

JS>1 

IF (THETA .GT. COSRAY) JS-KCRcST 

* 

♦* CALCULATE COEFFICIENTS OF CUBIC CURVE •* 
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160 Al-SLl 

A2»(3.0*DYPS-0XPS*(2.0*SL14SL2)»/0XPS**2 

A3»»-2.0*OYPS+OXPS*(SL1+SL2IJ/OXPS**3 

DX>DXPS/6.0 

NS-0 

* 

** CALCULATE AND DRAto POINTS ON CUBIC CURVE BETWEEN TWO RAYS ** 

* 

IF (JS .GE. 7) GO TO 180 

DO 170 J"4S»7 

NS-NS+1 

XPL(NS)«(XPL0T1 (!)♦( J-1)A0X)*AM/PLTSCAL 
D1-(J-1)*DX 

YPL(NS)-(YPL0T1(I)*01*(A1+D1*(A2+01*A3) n*AH/PLTSCAL 
IF (YPL(NS) .LT. 0.0) YPL(NS)-0.0 
170 CONTINUE 

CALL DRAW (XPL*YPL>NS) 

180 CONTINUE 

190 IF (NORAY .LT. 2) GO TO 210 

* 

♦♦ SAVE PREVIOUS TWO RAYS ♦* 

* 

DC 200 I>NFRST1>NLAST1 
XPL0T2m ■ XPLOTKI) 

200 YPL0T2(I) • YPLQTKI) 

210 DO 220 I-NFRST.NLAST 
XPLOTKI) ■ XPLOT(I) 

220 YPLOTKI) • YPLOT( I) 

NFRST2 - NFRSTl 
NLAST2 « NLASTl 
NFRSTl • NFRST 
NLASTl • NLAST 
GO TO 240 

* 

♦♦ END PLOTTING ♦♦ 

A 

230 CALL CALPLT(0.»0.»999) 

240 return 
END 



Subroutine RAYOPT 


This subroutine determines the initial ray spacing required to satisfy 
specified final (nearshore) ray spacing conditions. Final ray spacing is 
defined as the distance between the shorewardmost synchronous points on adja- 
cent rays. The routine would normally be used in conjunction with the wave- 
crest plotting option in order to facilitate the construction of curved crest 
segments between adjacent rays. A detailed discussion of this routine and 
pertinent input parameters is given in appendix C. The flow diagram and list- 
ing of subroutine RAYOPT follow. 
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RETURN 



55 





















SUBROUTINE RAYOPT (IREPEAT) 


^^^t^^^^4**t******************************************************^**** 

* THIS SUBROUTINE ITERATES FOR THE INITIAL RAY SPACING SO THAT 

* FINAL RAY SEPARATION CONDITIONS ARE SATISFIED . 

^t********************************************************************* 

COMMON Bl« 32» CO«CX Y# CPRIME« DCON> DELTAS# DEP(A2« 32 )#DISTMAX«D1STM IN* 
•^DRC#OTGR#DXY#GRID#GRINC#HO#IBLK(101)>IGO#IOUTUNT#IREF#ISTOP* JGO# 
^JPRTFRQ(3)#KCREST#KPL0T«LIMNPT#MI#MJ#NCREST#NFRST#NFRST1#NLAST# 
•»^NLAST1* NORA Y>NOROW#NPT>NSURFCE#NTIC#NTICaFF#NURITE# POX# PDY#PD2X> 
♦P02Y#PD2XY#PLTLNGX#PLTLNGY#PLTSCAL#RAYAN6#RAYMAX#RAYMIN#RAYSPC# 
>RCCO#RK#SK#SIG#TI0E#TITL(e)#V#WL#UL0#XO#XN#XPLOT(50O)#XPLOTl(5OO)# 

♦ XPL0T2(500),XTEST» Y0» YM, YPLOT ( 500 ) # Y PLOTl ( 500 ) # YPLOT2 { 500 )# YTEST 

IF (NORAY .LT. 2) GO TO BO 
NTEST ■ NLAST 

IF (NLAST .GT. NLASTl) NTEST-NLASTl 

* 

** COMPUTE FINAL SPACING BETWEEN ADJACENT RAYS ** 

* 

DIST • SORT((XPLOT(NTEST)-XPLOT1(NTEST))**2+(YPLOT(NTEST)- 

♦ YPLOTKNTEST) )**2I 

IF (MINDIST .GT. 0) GO TO 30 
IF (DIST .LT. DISTMAX) GO TO 20 

* 

** IF DISTANCE BETWEEN RAYS IS TOO LARGE# REDUCE INITIAL RAY ** 

** SEPARATION AND REPEAT. MAX NO OF PASSES EQUAL 10 

* 

MAXDIST • 1 

IF (NPASS .GT. 10) GO TO 80 

* 

** IF INITIAL RAY SPACING LESS THAN MINIMUM# SET TO MINIMUM ♦♦ 

* 

IF (RAYSPC .LE. RAYMIN) GO TO 10 

RAY2 • RAYSPC 

DIST2 • DIST 

RAYSPC - 0.5*RAYSPC 

GO TO 70 

10 RAYSPC • RAYMIN 
GO TO 80 
20 MINDIST ■ 1 
NPASS-1 

30 IF (NPASS .GT. 10) GO TO 80 

IF (DIST .LT. OISTMIN) GO TO -^0 
IF (DIST .LE. DISTMAX) GO TO 80 
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RAY2 ■ RAYSPC 
DIST2 - OIST 
MAXOIST > 1 
GO TO 50 

* 

** IP PINAL DISTANCE BETWEEN RAYS IS TOO SNALL* PINO INITIAL RAY 

** SEPARATION BY NEWTON-RAPHSON METHOD AND REPEAT. MAX NO 

** OP PASSES EQUAL 10 

* 

AO RAYl ■ RAYSPC 
DISTl • OIST 

IP (NAXDIST .GT. 0) GO TO SO 
RAYSPC • 2.0*RAYSPC 
GO TO 60 

50 DISTM ■ 0.5*{DISTMIN*0ISTMAX» 

RAYSPC ■ RAY1-K0ISTM-0IST1)*(RAY2-RAY1)/(DIST2-DIST1) 

* 

** IP INITIAL RAY SPACING GREATER THAN MAXIMUM* SET TO MAXIMUM ** 

♦ 

60 IP (RAYSPC .LT. RAYMAX) GO TO 70 
RAYSPC > RAYMAX 
GO TO 80 

70 XTEST • XTESTL 
YTEST ■ YTESTL 
NTIC • NTICL 
NPASS - NPASSAl 
IREPEAT • 1 
GO TO 90 

* 

** GOOD PASS* RESET AND UPDATE INITIAL CONDITIONS ** 

* 

BO MAXOIST • 0 
IREPEAT ■ 0 
MINOIST • 0 
XTESTL • XTEST 
YTESTL ■ YTEST 
NTICL ■ NTIC 
NPASS ■ 1 

90 RETURN 
END 



Subroutine MATRIX 


« 


This is an auxiliary routine incorporated to perform the matrix multipli- 
cation required in the determination of the coefficients for the topography 
approximating polynomial used in subroutine DEPTH. The subroutine multiplies 
the invariant matrix SX which depends only on the topography approximating 
technique selected, by a column vector of depth values D to obtain the 
resultant column vector of coefficients C to be used with the particular 
polynomial. The flow diagram and listing of subroutine MATRIX follow. 



SUBROUTINE MATRIX ( IR> JC* SX>0> C ) 

t********************** *************** ******************************** 

THIS SUBROUTINE MULTIPLIES A MATRIX TIMES A VECTOR. 

♦♦ INPUT ♦* 

IR - NO OF ROWS IN MATRIX A 

JC - NO OF COLUMNS IN MATRIX A 

SX(IR.JC) ■ INPUT MATRIX 
0(JC) • INPUT VECTOR 

CUR) ■ OUTPUT VECTOR 

********************************************************************** 

DIMENSION SX(IR» JC )*D(16).C(16) 

DO 10 K-l#16 
10 C(K) - 0. 

00 20 I-l.IR 
DO 20 J-l.JC 

20 C(I) • C(I)+SX(I,J >*0(J) 

RETURN 

END 
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PROGRAM USAGE 


The computer program is coded in FORTRAN Extended (FTN) Version 4 for use 
on the Control Data Cyber 173 or 175 computers under NOS 1.2. A central memory 
field length of 660008 locations is required, along with a peripheral mass- 
storage file for input of bathymetry data as discussed subsequently. Central 
processing unit (CPU) time requirements for program execution vary with the 
size of the geographic region being considered, the spatial density of wave 
rays being computed within the region, and the graphic output mode selected. 


Bathymetry Data Input 

As a first step, a bathymetry data array corresponding to a square grid 
with a uniform spacing of GRID feet per division must be developed for the 
geographic region to be investigated. It is recommended that the data array 
be developed by using a special transverse Mercator map projection centered 
within the region of interest in order to minimize error in lengthy ray paths 
due to curvature of the Earth. Advantages of using such a projection are dis- 
cussed in reference 6, and the projection itself is discussed in detail in ref- 
erence 9. For input to the computer program, the regional bathymetry array 
must be divided into overlapped modules of dimensions 42 by 32. Modularization 
must be performed as a separate function and can be done using the utility 
program CREMOD, which is described in detail in appendix D. During program 
execution, the modularized bathymetry data file TAPE1 is initialized as a 
random-access file in main program WAVE by using the system subroutine OPENMS. 
Modules are then read from the file to central memory as required in subroutine 
DEPTH by using the system subroutine READMS. 


Card Input 

Other parameters required for program execution are input in card form by 
using standard Control Data NAMELIST. The program is designed to compute a 
number of sets (NOSETS) of refraction diagrams with a single computer run. 
Parameters which are fixed for all sets of diagrams to be computed are input 
under the NAMELIST group name NPUT1 and are listed with their default values 
in table I. If the user desires to use default values, only the parameters 
specifying the modular bathymetry grid configuration must be supplied. Param- 
eters in effect for a single set are input under the NAMELIST group name NPUT2 
and are listed with their default values in table II. If the user desires to 
use default values with this group, only the wave period and initial wave 
height and direction must be specified. In addition, a title card of up to 
80 alphanumeric characters is required before each group of NPUT2 cards as 
identification of the conditions in effect for that set. Normal termination 
of program execution and closure of the plot file is effected after NOSETS of 
ray diagrams have been computed. However, in the event of user oversight, an 
end of file (EOF) is checked on card input on both NPUT1 and NPUT2 to insure 
proper termination. 
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Output Description 


At the beginning of the printed output, the basic program input parameters 
are identified and the values assigned to them are printed. Omission by the 
user of any mandatory input parameter causes the program to print an appropri- 
ate error message and then terminate. If KPLOT is specified as 2 or 3, the 
additional input parameters required for those plotting options are identified 
and their values printed. At the start of printed values for each ray, a 
header is printed which gives the values of constants valid for the current ray 
set along with title information for the current ray. Computed ray parameters 
are printed at specified increments in the ray computational point index NPT; 
these increments are input in the array JPRTFRQ. Output units for spatial 
coordinates (x,y) are in nautical miles. Output units for other computed 
parameters can be either U.S. Customary Units or SI Units, as specified by the 
input parameter lOUTUNT. The computed parameters are printed in column format 
in the following order: 

NPT ray computational point index 

X alongshore coordinate x, nautical miles 

Y offshore coordinate y, nautical miles 

A ray direction a, degrees 

DXY water depth d, feet or meters 

WL wavelength L, feet or meters 

C wave speed c, feet/second or meters/second 

RK refraction coefficient Kp 

SK shoaling coefficient Kg 

H wave height H, feet or meters 

At the end of printed parameters for each ray, the condition which led to ray 
termination is printed. 

The plotted output, depending upon the value of KPLOT, gives either ray 
diagrams or crest diagrams . The length of each axis is determined by the input 
parameters PLTLNGX and PLTLNGY. Scaling of the data and the axis scales are 
controlled by the parameter PLTSCAL. Major tick marks are added at 10-nautical 
mile intervals, and minor tick marks are added at 1 -nautical -mile intervals 
along both axes. In plotting ray diagrams, a plus symbol (+) is drawn at incre 
ments of NCREST computational points along each ray, and every fifth ray is 
numbered. Initial ray spacing is specified by RAYSPC in the default mode or 
is computed to meet constraints on final (nearshore) ray spacing if that option 
is selected. 


60 



Sample Cases 


In order to illustrate the input and output options of the program, three 
sample cases are shovm. The cases represent ray and crest diagrams for the 
Baltimore Canyon region of the mid-Atlantic continental shelf, extending from 
near Wachapreague Inlet, Virginia, to near Manasquan Inlet, New Jersey. The 
bathymetry data array for this region was divided into a total of 72 modules 
of dimensions 42 by 32, with 8 rows and 9 columns of data modules. For each 
sample case, the basic program parameters, the input parameters for the ray 
set, and the computed parameters for a representative ray within the set are 
shown. Also the ray or crest diagram constructed for the set in question is 
presented . It should be noted that the land masses and isobaths have been 
added to the basic diagrams for display purposes. 

Case 1 .- For this case, a ray diagram was generated with equal initial 
spacing between the computed rays. KPLOT was specified as 1, and lOUTUNT was 
specified as 1 , resulting in output parameters other than x and y being 
given in U.S. Customary Units. This case required 10 seconds for computer 
processing. Figure 3 gives the ray diagram for sample case 1. A listing of 
the input parcimeters for sample case 1 is given as follows: 

NPUT1 NPUT2 


DCON 

1.0 

DELTAS 

. . 0.01 

GRID 

. . 3038.06 

GRINC 

1.0 

lOUTUNT . . . . 

. . 1 

JPRTFRQ(I) . . . 

. . 20 

JPRTFRQ(2) . . . 

. . 10 

JPRTFRQ(3) . . . 

. . 5 

KPLOT 

. . 1 

LIMNPT 

. . 400 

MI 

. . 300 

MJ 

. . 254 

NOMOD 

. . 72 

NOROW 

. .• 8 

NOSETS 

. . 1 

PLTLNGX . . . . 

. . 12.0 

PLTLNGY . . . . 

. . 10.0 

PLTSCAL . . . . 

. . 12.5 

TIDE 

. . 0.0 


A 

. . . -45.0 

HO 

... 6.0 

NCREST 

... 5 

NSURFCE . . . . 

... 1 

RAYSPC 

... 5.0 

T 

. . . 12.0 
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The output for sample case 


1 follows: 


BASIC PROGRAM PARAMETERS 



DEPTH UNIT MULT... 

.DC ON 

■ 

1.00 

MIN STEP SIZE.... 

.DE LTAS 

■ 

.01 

GRID SIZE. 

.GR ID 

m 

3038.0b 

DEEP WATER STEP SIZE 

.GRINC 

m 

1.00 

OUTPUT UNIT OPTION 

.I3UTUNT 

m 

1 

PRINT FREQ 1 (D/L>0.5) 

.JPRTFRQ(l) 

m 

20 

PRINT FREQ 2 ( 0.25<D/L<0. 5 ) . . 

.JPRTFRO(2) 

m 

10 

PRINT FREQ 3 (D/L<0.25) 

.JPRTFRQU) 

m 

5 

PLOT OPTION.. 

.KPLOT 

m 

1 

MAX NO OF RAY POINTS 

.LIMNPT 

m 

400 

GRID LIMITS ABSCISSA 

.MI 

m 

300 

ORDINATE....... 

.MJ 

m 

254 

NO OF DEPTH MODULES 

• NOMOD 

m 

72 

NO OF ROW MODULES IN NOMOD... 

. NOROW 

m 

8 

NO OF SETS OF RAYS 

•NOSETS 

m 

1 

PLOT X AXES LENGTH 

.PLTLNGX 

m 

12.00 

PLOT Y AXES LENGTH 

•PLTLNGY 

m 

10.00 

PLOT SCALE FACTOR 

.PLTSCAL 

m 

12.50 

TIDE HEIGHT 

.TIDE 

m 

0.00 
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SAMPLE CASE NO 1 

PERIOD • 12.00 SEC. ..TIME STEP - ^9.44 SEC. ..RAY SPACING - 5.0. ..CREST SPACING - 5.. .BOTTOM APPROX CODE 

PAY NO - 2. ..SET NO ■ 1 


k 


X 


o 

o ^ 



o 

O' 

r4 

s> 

OJ <D ^ 

<Z> 

lO 

^ xj* tA 

o 

o 

m 



xT 


lA 

fO 

•4- 


O 

CM 

nO 

O 

lA 

lA 

<X> 


CO 

o 


o 

o « 


m 



o 


^ tn ^ 
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RAY STOPPED. REACHED SHORE 



Case 2.- This case illustrates a ray diagram in which the final (nearshore) 
ray separation distance is controlled. KPLOT was set equal to 2 and the limit- 
ing conditions on initial and final ray separation distances are printed at the 
end of the program parameter list. The output units are in SI Units as speci- 
fied by setting lOUTUNT equal to 2. This case required 22 seconds for computer 
processing. Figure 4 shows the ray diagram for this sample case. The listing 
of the input parameters for sample case 2 is given as follows: 

NPUT1 NPUT2 


DCON 

1 .0 

A 

. . . -45.0 

DELTAS 

0.01 

HO 

... 6.0 

GRID 

3038.06 

NCREST 

... 5 

GRINC 

1 .0 

NSURFCE . . . . 

... 1 

lOUTUNT 

2 

RAYSPC 

... 1.0 

JPRTFRQ(I) 

20 

T 

. . . 12.0 

JPRTFRQ(2) 

10 

DISTMAX . . . . 

... 50.0 

JPRTFRQO) 

5 

DISTMIN . . . . 

. . . 30.0 

KPLOT 

2 

KCREST 

... 1 

LIMNPT 

400 

RAYANG 

. . . 120.0 

MI 

300 

RAYMAX 

. . . 15.0 

MJ 

254 

• RAYMIN 

. . . 0.01 

NOMOD 

72 



NOROW 

8 



NOSETS 

1 



PLTLNGX 

12.0 



PLTLNGY 

10.0 



PLTSCAL 

12.5 



TIDE 

0.0 
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The output for sample case 2 follows: 


BASIC PROGRAM PARAMETERS 


DEPTH UNIT MULT 

MIN STEP SIZE............... 

..DC ON 
. .DELTAS 

• 

m 

1.0 

.01 

3038.0 

1.00 

2 

GRID SIZE.... 

..GRID 

m 

DEEP WATER STEP SIZE........ 

. .6R INC 

m 

OUTPUT UNIT OPTION.......... 

. . I3UTUNT 

m 

PRINT EREO 1 I D /I >0. 5 1 . . . . . . 

• . JP RTFRQI 1 ) 

m 

20 

10 

PRINT FREQ 2 <0.25<0/L<0. 5). 

..JPRTFRQ(2) 

m 

PRINT FREQ 3 (0/L<0* 25) 

. . jp RTFR0(3 ) 

m 

5 

PLOT OPTION................. 

. .KPLOT 

m 

2 

MAX NO OF RAY POINTS........ 

. . LIHNPT 

m 

<>00 

300 

25A 

72 

8 

GRID LIMITS ARSCISSA. ..... 

. .Ml 

m 

ORDINATE. .... . 

• .MJ 

m 

NO OF DEPTH MODUt FS ........ . 

• .NO MOD 

m 

NO OF ROW MODULES IN NOMOO*. 

..NOROW 

m 

NO OF SETS OF RAYS.......... 

. .NOSETS 

9 

1 

Pi OT X AXES 1 FNGTH. ........ . 

. .PLTLNGX 

m 

12.00 

10.00 

12.50 

0.00 

50.0 

30.0 
1 

PLOT Y AXES LENGTH.......... 

. .PLTLNGY 

m 

PLOT SCALE FACTOR........... 

. .PLTSCAL 

m 

TIDE HEIGHT................. 

. .TI DE 

9 

MAX FINAL RAY SPACING....... 

. .DISTMAX 

m 

MIN FINAL RAY S PAC I NG . . . . . . . 

. .DI STMIN 

9 

INITIAL PT FOR CREST CURVE.. 

..KCREST 

9 

MIN ANGLE BETWEEN SEGMENTS-.. 

..RAYANG 

9 

120.0 

MAX INITIAL RAY SPACING.^jjj 

• .RAYMAX 

9 

15.0 

.01 

MIN INITIAL RAY SPACING 

..RAYMIN 

9 
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SAMPLE CASE NO 2 

PERIOD - 12.00 SEC...rr«£ STEP - 49.44 SEC. ..RAT SPACISG • 9,0,,, CREST SPACING ■ .BOTTOM APPROX CODE 
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Case 3 .- The crest diagrams constructed from the rays computed in sample 
case 2 are illustrated by this sample case. KPLOT was set to 3 to specify this 
plotting option. The output units are U.S. Customary as specified by setting 
lOUTUNT to 1. The resulting crest diagram is shown in figure 5. This case 
required 22 seconds for computer processing. The listing of the input param- 
eters for sample case 3 is given as follows: 

NPUT 1 NPUT2 


DCON 

1 .0 

A 

. . . -45.0 

DELTAS 

. . 0.01 

HO 

... 6.0 

GRID 

. . 3038.06 

NCREST .... 

... 5 

GRINC 

1.0 

NSURFCE . . . 

... 1 

lOUTUNT . . . . 

. . 1 

RAYSPC .... 

... 1.0 

JPRTFRQ(I) . . . 

. . 20 

T 

. . . 12.0 

JPRTFRQ(2) . . . 

. . 10 

DISTMAX . . . 

. . . 50.0 

JPRTFRQ(3) . . . 

. . 5 

DISTMIN . . . 

. . . 30.0 

KPLOT 

. . 3 

KCREST .... 

... 3 

LIMNPT 

. . 400 

RAYANG .... 

. . . 120.0 

MI 

. . 300 

RAYMAX .... 

. . . 15.0 

MJ 

. . 254 

RAYMIN .... 

. . . 0.01 

NOMOD 

. . 72 



NOROW 

. . 8 



NOSETS 

. . 1 



PLTLNGX . . . . 

. . 12.0 



PLTLNGY . . . . 

. . 10.0 



PLTSCAL . . . . 

. . 12.5 



TIDE 

. . 0.0 
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The output for sample case 3 follows: 


BASIC PROGRAM PARAMETERS 


DEPTH UNIT MULT*. 

MIN STEP SIZE.**.*...******... 

GRID SIZE 

DEEP WATER STEP SIZE 

OUTPUT UNIT OPTION 

PRINT FREQ 1 ( D / L>0« 5 . 

PRINT FREQ Z (0*Z5<0/L<0. S) . . • 
PRINT FREQ 3 (D /L<0*25) ...... . 

PLOT OPTION*. 

MAX NO OF RAY POINTS ****••••• • 

GRID LIMITS ABSCISSA 

ORDINATE... 

NO OF DEPTH MODULES ***•*...*. * 
NO OF ROW MODULES IN NOMOD.... 

NO OF SETS OF RAYS 

PLOT X AXES LENGTH. .*•*.• 

PLOT Y AXES LENGTH... 

PLOT SCALE FACTOR 

TIDE HEIGHT. 

MAX FINAL RAY SPACING**** **•• * 

MIN FINAL RAY SPACING 

INITIAL PT FOR CREST CURVE.... 
MIN ANGLE BETWEEN SEGMENTS.... 

MAX INITIAL RAY SPACING 

MIN INITIAL RAY SPACING 
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APPENDIX A 


DESCRIPTION OF SYSTEM SUBROUTINES 

The program described in this paper uses several system subroutines to 
satisfy various data handling and plotting requirements. This appendix lists 
the names of subroutines and their functions. 

AXES Draws a line, annotates the value of the variable at specified 

intervals with or without tick marks, and provides an axis 
identification label. 

CALPLT Creates a reference point for plotting and ends plotting. 

CLOSMS Closes the random-access mass-storage file. 

DRAW Draws a line between specified points. 

ENCODE Reformats data in memory. 

LEROY Sets plot speed . 

NFRAME Advances the plot frame and sets the plot origin. 

NOTATE Draws alphanumeric information for plot annotation and labeling. 

NUMBER Converts a floating point number to binary coded decimal (BCD) 

format and draws the resulting alphanumeric characters. 

OPENMS Initializes the random-access mass-storage file. 

PNTPLT Draws a given symbol centered on a given coordinate value. 

PSEUDO Initializes and writes a plot vector file. 

READMS Transfers data from mass storage to central memory. 

WRITMS Transfers data from central memory to mass storage. 
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CONSTRUCTION OF WAVE-CREST DIAGRAMS 

The interaction of ocean waves with coastal bottom topography can be 
examined most readily by plotting wave-crest patterns rather than ray patterns. 
Wave-crest patterns can be constructed by connecting the same point in time on 
adjacent rays by a curve which is locally orthogonal to each of the rays it 
intersects. The simplest smooth curve which can be drawn perpendicular to two 
adjacent rays is represented by a third-order equation. A third-order curve 
passing through point j on ray i-1 can be written 


y = yi-1,j -*■ - ^i-1,j) + ^2(x - + 33(x - Xi_i j)3 (B1) 

where a-| , a 2 , and a 3 are coefficients to be determined. The coefficients 
can be found by constraining the curve to pass also through point j on ray i 
and constraining the slope of the cubic curve (dy/dx) to be equal to the normals 
at point j on rays i-1 and i (fig. 6(a)). The slopes required for the 
curve to be orthogonal to rays i-1 and i at point j can be approximated 
by 


/dy\ ^i-1,j-1 ■ ^i-1,j+1 

\dx/i_i J yi- 1 ,j +1 - yi- 1 ,j -1 


(B2) 


/dy\ ^i,j-1 ■ ^i,j+1 

Vdx/ij yi,j+1 - yi,j-1 


With these results, the three coefficients in equation (B1) are determined to 
be 


ai 



(B4) 



(B5) 
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^3 


. (5:') 


(xi,j - xi_ij)3 




(B6) 


The cubic equation is used to determine the points on the crest segment between 
two adjacent rays. To give the crest curve a smooth appearance, seven points 
are used to describe the crest segment. 

When waves approach intermediate water and refraction begins to occur, 
adjacent rays can cross. Straightforward construction of crest segments in 
such situations is probably unrealistic in a physical sense and could often 
result in a zigzag crest pattern. To circumvent this problem, logic has been 
added to delete all or part of the crest segment between crossing adjacent 
rays i-1 and i. In order to signal the onset of this situation, the angle 
between adjoining first-order crest segments is computed. In figure 6(b), the 
angle 6 between the first-order segment joining point j on ray i-2 to 

point j on ray i-1 , and the first-order segment joining point j on 

ray i-1 to point j on ray i is given by 


0 = cos“^ 


(Xj - Xj.i)(xi.2 - -h (yj - yj-i)(yi.2 - Yj-i) 

VUi - Xi _-()2 + (y^ - yi_i)2 \/(Xi_2 - 


(B7) 


Whenever this angle is less than a specified input parameter (RAYANG) , the 
curved crest segment connecting point j on ray i-1 to point j on ray i 
is modified. Since a maximum of seven points are used to describe the crest 
segment, the modification is accomplished by starting the segment at an inter- 
mediate point between the rays rather than at point j on ray i-1 . This 
results in crest patterns being disjointed when adjacent rays i-1 and i 
cross. The program input parameter specifying the point at which to begin con- 
struction of the modified curved crest segment is KCREST and can be any number 
between 1 and 7 inclusive. The disjointed effect is illustrated in figure 5 
with RAYANG equal to 120° and KCREST equal to 3 • In this case only two-thirds 
of the crest curve is drawn when the angle between adjoining first-order crest 
segments is less than 120°. The logic for this procedure is incorporated in 
subroutine PLOTR. Crest plotting is initiated when the input parameter KPLOT 
is equal to 3 . 
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DISCUSSION OF SUBROUTINE RAYOPT 

When wave rays approach shallow water and refraction occurs, cases of 
extreme convergence or divergence can occur, resulting in either large areas 
completely devoid of rays or small areas covered by many rays. To allow the 
user a degree of control over nearshore ray spacing, the subroutine RAYOPT was 
developed. In subroutine RAYOPT, desired nearshore ray spacing is achieved by 
perturbation of the initial ray spacing. Upper and lower bounds specified by 
input parameters are imposed on both the initial ray spacing and the final ray 
spacing. If desired results are not obtained in 10 iterations, then the 10th 
iteration is specified as the solution. 

To start the RAYOPT procedure, an initially small value is input or 
assigned by default to the input parameter RAYSPC, which is the initial ray 
spacing. The initial ray in the set is found by using subroutine RAYINIT. A 
second ray is begun at a distance RAYSPC from the initial ray, its complete 
path is generated, and its final separation distance with respect to the pre- 
ceding ray is computed. If this distance is greater than the input parameter 
DISTMAX, the current RAYSPC is halved and a new trial ray is generated. Con- 
versely, if this distance is less than the input parameter DISTMIN, the current 
RAYSPC is doubled and a new trial ray is generated. This procedure continues 
until the final separation distance lies between the limits DISTMIN and DISTMAX. 
If at anytime in the iteration process RAYSPC falls below the input value RAYMIN 
or exceeds the input value RAYMAX, the iteration process terminates and RAYSPC 
is set to the respective limit. In some cases, the final ray separation dis- 
tance oscillates between values of DISTMIN and DISTMAX. When this occurs, an 
acceptable initial ray spacing is computed by using a Newton-Raphson iteration 
scheme. For all cases, the final iterated value for RAYSPC is used as the 
first trial value in computing the next ray. 

Graphical results obtained by using subroutine RAYOPT can be seen by com- 
paring figures 3 and 4, Figure 3 presents a ray diagram computed with a uniform 
initial ray spacing (RAYSPC) of 5 grid units, with no control exercised over 
final spacing. There are large gaps between adjacent rays in the nearshore 
region in this diagram. In figure 4, final ray spacing is controlled with 
DISTMIN equal to 30 grid units and DISTMAX equal to 50 grid units. The initial 
ray spacing is constrained to lie between 0.01 grid units (RAYMIN) and 15 grid 
units (RAYMAX). As can be seen, the nearshore ray spacing is much more uniform 
than that in figure 3. In deeper areas in which no refraction occurs, ray 
spacing is increased by the subroutine until either RAYMAX or DISTMAX is 
exceeded. User control over the final ray spacing is particularly useful in 
cases in which crest curves are to be constructed between adjacent rays, as 
described in appendix B. Figure 7 shows the complete crest pattern constructed 
for the ray diagram shown in figure 3, which was generated by using a constant 
initial ray spacing of 5 grid units. In constructing the crest patterns of 
figure 7, RAYANG was specified as 120° and KCREST was specified as 1. This 
crest pattern does not accurately depict nearshore wave fronts because of the 


73 



APPENDIX C 


sparsity of wave rays in certain nearshore areas. The crest pattern of fig- 
ure 5 was constructed for the ray diagram of figure 4, in which final ray spac- 
ing was controlled. Here, a more realistic version of the wave fronts for the 
nearshore zone is seen. 


An experimental period may be required on the part of the user to obtain 
an acceptable crest pattern for a particular application. Care should be taken 
when selecting values for the input parameters required in generating crest 
patterns. The following discussion on each of the required input parameters 
is included as a guide in specifying these parameters for any application: 


DISTMAX 


DISTMIN 


KCREST 


RAYANG 


RAYMAX 


RAYMIN 


This input constant is the maximum permissible final ray separa- 
tion distance. Too large a number causes wide nearshore ray 
separation, such as that seen in figure 3, and results in 
unsatisfactory crest patterns (fig. 7). Too small a number 
causes too many rays to be generated and results in patterns 
similar to that shown in figure 8, where DISTMAX is equal to 
30 grid units. In order for the iteration process to converge 
in a few cycles, this number should differ from DISTMIN by a 
reasonable margin. 

This parameter is the minimum permissible final ray separation 
distance. This number should be large enough to insure suffi- 
cient nearshore separation between adjacent rays. 

This parameter controls modification of the curved crest segment 
to be drawn between adjacent rays which cross. Normally, a 
value equal to 2 or 3 produces an acceptable diagram. Failure 
to draw at least part of the curved segment (i.e., setting 
KCREST =7) can result in large gaps in the pattern, as seen 
in figure 9 (compared with fig. 5). 

This parameter is the limiting angle between adjoining first- 
order crest segments. If, in constructing a crest pattern, the 
angle between adjoining first-order crest segments is less than 
RAYANG, construction of the current curved segment begins at 
point KCREST rather than at the first point. A value for RAYANG 
around 120° should give acceptable patterns. 

This parameter is the maximum permissible initial ray spacing. 
RAYMAX should be set to a large value to allow an optimum search 
for an acceptable initial ray point. In addition, only a few 
rays are required to define the crest pattern in situations in 
which the rays remain in deep water over their entire path. 

This parameter is the minimum permissible initial ray spacing. 

In all instances, this number should be set to an arbitrarily 
small value in order that final ray constraints can be met. 
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UTILITY PROGRAM CREMOD 

In the wave refraction computer model described in this paper, bathymetry 
data are stored and retrieved in the form of overlapped modules of dimensions 
42 by 32 . Random-access techniques are used to pass data modules as required 
from a peripheral mass-storage file to the central memory of the computer. The 
potential user, however, generally has bathymetry data available in the form 
of an array of dimensions MI by MJ covering the entire region of interest. 

This appendix describes a utility program CREMOD which can be used to modu- 
larize the regional array supplied by the user. Upon completion of execution 
of this program, the modularized data reside on the peripheral file TAPE1, 
which then can be saved as a permanent file by the user for later input to the 
refraction model. 

It is assumed that the regional array is supplied in card deck form and 
that the data are arranged in columnwise order with the first column corre- 
sponding to the x-axis of the region. The total number of rows of data MI and 
the total number of columns of data MJ are required input parameters. The user 
must also designate the format (FORMT) with which the data cards are read. To 
avoid possible extraneous results, it is strongly recommended that the format 
be chosen so that the input record length is an integral divisor of MI. The 
regional array is divided into a total of NOMOD overlapped modules which are 
indexed by columns, as shown in figure 10, and written to the disk file TAPE1 
by using system subroutine WRITMS. If the regional configuration is such that 
portions of certain modules extend beyond the boundaries of the MI by MJ array 
(as illustrated in fig. 10), data values within those portions are set equal 
to the values in the corresponding column of row MI or in the corresponding row 
of column MJ, as the case requires. Values for NOROW, the number of module 
rows, NOCOL, the number of module columns, and NOMOD are printed when modulari- 
zation is complete . 

Program CREMOD uses the system subroutines OPENMS, WRITMS, and CLOSMS for 
initializing, writing records to, and closing the mass-storage file TAPE1. Ref- 
erence 10 provides a detailed description of each of these routines. On the 
Control Data Cyber 175 computer under NOS 1.2, program CREMOD required a central 
memory field length of 640000 and, for a sample case in which MI and MJ were 
both 100, required 4 seconds of CPU time. The user should again be reminded 
to save the file TAPE1 as a permanent file after execution of program CREMOD. 

The flow diagram and listing of program CREMOD follow. 
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PROGRAH CREMOD ( INPUT> OUTPUT* TAPE 5«INPUT# TAPE 1 ) 

OIHENSION DEP(A20* 3 2) * OHOOC A2> 32) * 19LK ( 101 ) 

* 

READ (5#180) HI*HJ>FORHT 

* 

* DETERMINE NO. OF MODULE ROUS AND COLUMNS. 

* 

NOCOL ■ MJ/29 
NCOUT - MOO (MJ»29) 

IF <NCOUT .LT. 4) GO TO 10 
NOCOL • NOCOL^l 
GO TO 20 

10 NCOUT • NCQUT+29 
20 NOROW ■ MI/39 

NROUT • MOO (MI >39) 

IF (NROUT .LT. <.) GO TO 30 
NOROW ■ N0R0W>1 
GO TO 90 

30 NROUT - NR0UT«39 
90 NOMOO • NOROWPNOCOL 

A 

A OPEN MASS STORAGE FILE> TAPEl. 

A 

CALL OPENMS ( 1> IBLK> NOMOOaI.O ) 

A 

A READ 3 COMPLETE COLUMNS OF DATA. 

A 

READ (5>F0RMT) ( (0 EP( I > J ) > I>1> HI ) > J>1 > 3 ) 

00 170 NC*1*NQC0L 

IF (NC .LT. NUCOL) 60 TO 70 

A 

♦ FINAL MODULE COLUMN. ONLY (NCOUT-3) SUB-COLUMNS OF DATA LEFT. 

A 

IF (NCOUT .EO. 32) GO TO 70 

READ (5/FORMT) ( (0 EP ( I* J ) * I«l> MI ) * J>9*NC0UT ) 

A 

A SET VALUES IN EXTERIOR COLUMNS - TO VALUES IN LAST COLUMN. 

A 

J1 • NCOUT+1 

00 60 J-J1.32 

00 50 I«l,MI 

DEP(I.J) • DEPd. NCOUT) 

50 CONTINUE 
60 CONTINUE 
60 TO 80' 

A 

A FULL MODULE COLUMN AVAILABLE. READ COLUMNS 9 THRU 32. 
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70 READ (5#F0RMT) ( (OEP ( I* J ) # I»l> MI )» J- A# 32) 

* 

* FILL DMOD ARRAY FOR SOROW MODULES IN THIS COLUMN, WRITE EACH TO 

* MASS STORAGE FILE. 

* 

80 DO lAO NR«l,NQRaW 
" IZERO ■ 39*(NR-1) 

NM ■ (NC-1)*N0R0W>NR 
IF (NR .LT. NOROW) GO TO 90 

* 

* FINAL MODULE ROW. NROUT-3 MORE SUB-ROWS OF DATA AVAILABLE. 

* 

IMAX - NROUT 
GO TO 100 
90 IMAX - A2 
100 DO 130 J>1,32 
DO 110 I-l.IMAX 
DM00(1,J) ■ OEP(IZERO«I,J) 

110 CONTINUE 

IF (IMAX .EQ. 42) GO TO 130 
11 • IMAX+1 
DO 120 1-11,42 
OMOOd.J) • DMOD(IMAX.J) 

120 CONTINUE 
130 CONTINUE 

CALL WRITMS ( 1, DMOD, 1344, NM) 

140 CONTINUE 

IF (NC .EO. NOCOL) GO TO 170 

« 

♦ SUB-COLUMNS 30-32 WILL BE SUB-COLUMNS 1-3 OF NEXT MODULE. 

* 

DO 160 J-1,3 
00 ISO I-1,MI 
OEP(I,J) • 0EP(I,J*29) 

150 CONTINUE 
160 CONTINUE 
170 CONTINUE 

PRINT 190, NOROW, NOCOL, NOMOD 
CALL CLOSMS (1) 

RETURN 

* 

180 FORMAT (215, AlO) 

190 FORMAT (2X,ANOROW I 3/2X, *NOCOL ■♦, I 3/2X, *NOMOD -*,13) 

END 
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TABLE I.- DESCRIPTION OF ELEMENTS IN NAMELIST NPUT1 


Variable 

name 

Description 

Default 

value 

DCON 

Multiplier to convert depth units to feet for 
internal computations 

1 .0 

DELTAS 

Minimum step length along ray in shallow water, 
grid units 

0.01 

GRID 

Number of feet per grid division 

None 

GRINC 

Step length along ray in deep water, grid units 

1 .0 

lOUTUNT 

Output units: nautical miles for spatial 

coordinates x and y; for all other parameters; 
1 = U.S. Customary Units; 2 = SI Units 

2 

JPRTFRQd ) 

Print frequency of ray parameters when depth is 
greater than one-half the wavelength, time steps 

20 

JPRTFRQ(2) 

Print frequency when depth is greater than one- 
fourth the wavelength but less than one-half the 
wavelength, time steps 

10 

JPRTFRQO) 

Print frequency when depth is less than one- 
fourth the wavelength, time steps 

5 

LIMNPT 

Maximum number of computation points along a ray 

400 

KPLOT 

Plotting option: 

0 = No plot 

1 = Plot wave rays with uniform initial spacing 

2 = Plot wave rays with controlled final spacing 

3 = Plot wave crests 

0 

MI 

Total number of rows in regional bathymetry data 
array 

None 

MJ 

Total number of columns in regional bathymetry 
data array 

None 

NOMOD 

Total number of 42 by 32 bathymetry data modules; 
maximum allowed is 100 

None 

NOROW 

Number of modules in x-direction 

None 

NOSETS 

Number of sets of rays to be computed 

1 

PLTLNGX 

Length of x-axis for plot, inches 

10.0 

PLTLNGY 

Length of y-axis for plot. Inches 

10.0 

PLTSCAL 

Scale factor for x- and y-axes; change in x or 
y value per inch of plot (input required for 
plotting only) 

None 

TIDE 

Height of tide, feet 

0.0 
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TABLE II.- DESCRIPTION OF ELMENTS IN NAMELIST NPUT2 


Variable 

name 

Description 

Default 

value 

A 

Initial ray direction, degrees 

None 

HO 

Deepwater wave height, feet 

None 

NCREST 

Spacing of crest tick marks or crest curves; 
integer multiple of program computational time 
step 

5 

NSURFCE 

Bottom topography approximation technique: 

1 = Quadratic least squares 

2 = Cubic least squares 

3 = Constrained bicubic interpolation 

1 

RAYSPC 

Initial spacing between wave rays in deep water, 
grid units 

5.0 

T 

Wave period, seconds 

None 

DISTMAXa 

Maximum permissible final separation distance 
between adjacent rays, grid units 

999. 

DISTMIN® 

Minimum permissible final separation distance 
between adjacent rays, grid units 

25. 

KCREST^ 

Point (between 1 and 7) between adjacent rays at 
which to begin plotting curved crest segment 
when angle between adjoining first-order segments 
is less than RAYANG 

1 

RAYANOa 

Minimum acceptable angle between adjoining first- 
order crest segments, below which crest plotting 
is modified, degrees 

120. 

RAYMAXa 

Maximum permissible initial ray spacing, grid 
units 

25. 

RAYMINa 

Minimum permissible initial ray spacing grid 
units 

1 

0.01 


^These parameters are used only for ray sets in which nearshore 
ray spacing is controlled (KPLOT = 2 or 3). 





Figure 1 Determination of starting point for first ray for ray direction 

greater than -90°. 
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